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L. _ STATEMENT  OF  WORK 

Conduct  a  theoretical  research  program  to  develop  quantum 
mechanical  methods  of  studying  nonadiabatic  effects  in  three-dimensional 
atom-diatom  collisions. 

II.  DESCRIPTION  OF  PROBLEM 


Chemical  dynamics  has  reached  the  stage  of  development  that 

allows  the  first-principles  determination  of  detailed  state-to-state 

information  for  many  kinetic  processes.^  These  new  experimental  and 

theoretical  methods  are  beginning  to  provide  information  of  significant 

importance  to  military  technology  in  such  areas  as  chemical  and  excimer 

laser  development,  studies  of  the  interaction  of  modern  weapons  systems 

with  the  atmosphere,  the  characterization  of  the  radiation  from  rocket 

(2  31 

plumes,  and  combustion  and  propellant  research.  ’  All  of  these  areas 
require  rate  data  for  specific  quantum  transitions  as  input  to  sophisticated 
kinetic  codes. 


2 


In  many  instances,  the  transitions  of  interest  are  difficult 

to  investigate  experimentally  due  to  short  lifetimes,  low  intensities, 

competing  processes,  or  simply  economic  factors.  Theoretical  approaches 

can  not  only  provide  information  of  extreme  utility  to  the  experimentalist 

as  support  for  the  interpretation  of  data,  but  in  their  own  right  can  be  the 

(41 

most  cost-effective  means  of  obtaining  such  information.  ' 

The  goal  of  this  research  program  is  to  develop  new  and  more 
efficient  quantum  scattering  methods  that  will  be  useful  in  applications 
to  state- to-state  collision  processes  involving  two-  and  three-atom  systems. 
The  emphasis  is  on  nonadiabatic  processes,  particularly  those  that  involve 
the  transfer  of  electronic  energy.  This  approach  is  based  on  the  coupled- 
channel  method,  and  stresses  reliable  approximations  that  allow  the  study  of 
light,  first-row  molecular  systems  involving  up  to  three  atoms.  An  important 
constituent  of  the  present  approach  is  the  incorporation  of 
potential  energy  surfaces  and  couplings  obtained  from  ab-initio  quantum 
chemistry.  The  suitable  analytical  representation  of  such  surfaces  is 
an  important  component  of  this  research. 

The  Born-Oppenheimer  (BO)  separation  of  electronic  and  nuclear 
motion  is  a  valuable  tool  in  molecular  theory  since  many  low-energy 
collision  processes  are  often  adequately  described  by  considering  motion 
on  a  single  potential  energy  surface.  For  inelastic  collisions  where 
avoided  crossings  or  small  separations  between  electronic  states  occur, 
and  for  reactions  which  involve  the  breaking  of  chemical  bonds  and 
reorganization  of  spin  couplings,  the  BO  approximation  can  be  a  poor 
one.  In  such  cases  it  is  necessary  to  consider  the  mixing  of  two  or  more 
adiabatic  electronic  states  that  arise  due  to  nuclear  motion. 
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Quantum  chemistry  is  concerned  with  solution  of  the  Schrodinger 
equation  for  electronic  motion  that  results  from  application  of  the  Born- 
Oppenheimer  approximation.  The  set  of  solutions  (adiabatic)  to  the 
electronic  problem  can  be  used  as  a  basis  for  expanding  the  total  wavefunctions 
for  nuclear  motion.  When  the  BO  approximation  is  valid,  the  wavefunctions 
for  nuclear  motion  adequately  describe  molecular  collisions  on  the  appropriate 
potential  surface  (neglecting  spin  effects). 

These  adiabatic  electronic  functions  can  often  be  strongly  coupled 

(51 

by  operators  neglected  in  the  BO  separation.  '  For  body-fixed  coordinates 
(BF),  in  which  the  electronic  problem  is  conveniently  solved,  these  operators 
take  the  form  of  BO  couplings  for  internal  motion  and  coriolis  couplings 
resulting  from  tumbling  of  the  BF  axis.  Spin-orbit  interaction  is  normally 
neglected  in  solving  for  the  electronic  eigenfunctions,  but  it  must  be 
included  with  the  nonadiabatic  couplings  for  a  proper  treatment  of  the  collision 
problem.  Other  terms  in  the  Breit-Pauli  hamiltonian  resulting  from  removing 
the  center-of-mass  motion  may  be  neglected  in  problems  of  chemical  interest. 

Since  most  quantum  coupled-channel  methods  make  use  of  partial - 

wave  expansions,  studies  of  interacting  open-shell  species  must  explicitly 

consider  the  various  angular  momentum  couplings  that  occur.  Several  quantum 

treatments  of  multiple-surface  effects  in  F  +  have  done  so.  Miller  and 

Wyatt^  and  DeVries  and  George^  utilize  the  valence  bond  character  of 

(8) 

DIM  theory  in  their  formulations,  while  Rebentrost  and  Lesterv  '  employ 
SCF  wavefunctions.  Depending  on  the  spin  and  angular  momentum  of  the 
collision  partners,  different  coupling  schemes  are  required.  These  studies 
are  the  only  ones  reported  for  the  interaction  of  a  structured  atom  with  a 
molecule  in  a  state.  Extensions  to  open-shell  molecules  are  necessary  to 
reach  the  ultimate  goal  of  this  project. 
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All  three  studies  mentioned  above  employ  diabatic  representations 

(51 

for  solving  the  coupled  equations.  '  These  are  usually  obtained  by  various 

prescriptions  from  the  adiabatic  representation,  and  are  not  unique.  They 

can  be  obtained  by  a  unitary  transformation  that  globally  eliminates  certain 

couplings.  The  advantage  of  diabatic  representations  is  that  one  can 

minimize  or  eliminate  the  large  couplings  due  to  nuclear  motion  and  instead 

employ  a  nondiagonal  representation  of  the  electronic  hamiltonian.  By 

eliminating  the  first  derivative  term,  the  coupled  equations  can  be 

integrated  using  very  efficient  algorithms.  Similar  couplings  appear  for 

vibrational  and  rotational  motions  in  reaction  coordinate  formulations  of 

reactive  scattering.  Although  the  equations  can  be  integrated  with  such 

(91 

terms  included,  better  stability  is  obtained  if  they  are  eliminated.  ' 

Since  ab  initio  adiabatic  potential  surfaces  and  their  couplings  are 
employed  in  this  approach,  it  is  important  that  adequate  methods  be 
developed  for  integrating  the  appropriate  coupled  equations. 

In  summary,  this  research  program  attempts  to  bring  together  the 
computational  tools  necessary  to  determine  from  first-principles,  state-to- 
state  probabilities  for  quantum  transitions  involving  rotational,  vibrational, 
and  electronic  degrees  of  freedom  for  atom-diatomic  molecule  collisions. 

In  Part  III  we  highlight  the  principal  objectives  of  this  program,  and  in 
Part  IV  we  examine  the  goals  achieved  during  the  past  18  months.  Part  V 
provides  recommendations  for  further  work. 
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III.  RESEARCH  OBJECTIVES 


The  overall  objectives  of  this  research  program  are  as  follows: 


•  Extend  the  3-D  reaction  coordinate  theory  of  chemical 
reactions  to  include  nonadiabatic  electronic  transitions. 

•  Investigate  various  decoupling  approximations  for 
reducing  the  complexity  of  the  coupled  equations. 

•  Develop  systematic  approaches  for  the  analytical 
representation  of  ab  initio  potential -energy  surfaces 
and  couplings. 

•  Develop  efficient  algorithms  for  the  integration  of 
coupled  equations  involving  nonadiabatic  couplings 
between  rotational,  vibrational,  and  electronic  states. 

•  Implement  these  methods  into  efficient  scattering  codes. 

•  Apply  these  codes  to  a  variety  of  problems  of  current 
interest. 


IV.  RESEARCH  ACCOMPLISHMENTS 


This  research  contract  was  originally  funded  as  a  36-month  effort, 
but  this  was  reduced  to  18  months  because  the  principal  investigator  left 
Battelle.  In  spite  of  the  short  duration  of  this  project,  there  have  been 
a  number  of  accomplishments  that  will  form  the  basis  of  a  practical  method 
for  studying  electronic  excitation  in  molecular  collisions. 

As  a  result  of  early  work  on  this  project,  rate  constants  from 
3-D  reactive  scattering  calculations  for  F  +  H£  and  H  +  H2 ( V=1 )  are  shown 
to  be  in  general  agreement  with  experiment  (see  Appendix  A).  Recent 
experiments  on  F  +  suggest  the  existence  of  a  resonance^)  predicted 
earlier  by  our  theoretical  approach. Surface  fitting  procedures  have 
been  developed  for  fitting  ab  initio  potential  energy  surfaces  (see 
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Appendix  B).  A  variety  of  techniques  for  integrating  coupled  equations 
were  investigated,  in  part  in  collaboration  with  the  NRCC  workshop  on 
computational  algorithms  in  scattering  theory.  These  methods  were  applied 
to  a  variety  of  vibrationally  and  rotationally  nonadiabatic  processes 
(Appendix  A)  and  to  electronically  nonadiabatic  processes  in  K  +  H 
collisions  (Appendix  C). 

Specific  accomplishments  are  as  follows: 

•  Previously  computed  reaction  probabilities  for  F  +  Ho 
were  used  to  determine  cross  sections  and  rate  constants 
for  this  reaction.  This  is  the  first  3-D  quantum 
mechanical  calculation  of  the  rate  of  a  chemical  reaction 
other  than  H  +  H2-  Arrhenius  parameters  from  the 
theoretical  calculations  are  in  reasonable  agreement  with 
experiment.  Perhaps  the  most  important  result  is  that 
it  is  possible  to  compute  probabilities  at  enough  values 
of  energy  and  total  angular  momentum  to  obtain  total 
state-to-state  cross  sections  over  the  range  of  energies 
required  to  compute  a  thermal  rate.  This  is  further 
discussed  in  Appendix  A. 

0  Quanta!  rate  constants  calculated  at  300  K  for  H  +  ^(IM) 
agree  with  some  experimental  results  and  are  in  apparent 
disagreement  with  classical  mechanics.  The  potential 
surface  used  in  this  study  is  not  very  reliable,  but  gives 
rates  for  V=0  in  good  agreement  with  experiment  and  other 
theoretical  values  (see  Appendix  A).  This  work  is  discussed 


I 
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in  some  detail  in  a  recent  review  by  Schatz.  Again, 
it  is  significant  that  these  calculations  are  possible 
with  modest  computing  resources.  The  H  +  H2(V=1) 
calculations  were  done  on  a  VAX  minicomputer! 

•  Methods  have  been  developed,  based  on  the  many-body 
approach  of  Murrell ,  to  obtain  analytical  representa¬ 
tions  of  three-  and  four-atom  potential  energy  surfaces. 
To  date  applications  have  been  made  to  0(^D)  +  H2(^2g), 
C(3P)  +  02(V),  0{3P)  +  H20  and  0(3P)  +  C02-  Quartic 
force  fields  for  H20  and  C02  are  accurately  reproduced 
with  this  technique  (see  Appendix  8). 

•  Codes  for  generating  3-D  electronic  correlation  diagrams 
in  reaction  coordinates,  including  rotational -vibrational 
degrees  of  freedom,  have  now  been  developed.  These  are 
general  codes  capable  of  treating  one  potential  surface 
at  a  time  and  are  not  restricted  to  linear  reaction 
intermediates.  An  analytical  representation  of  the 
potential  surfaces  for  each  electronic  state  is  required 
input.  Systems  studied  so  far  are  F  +  H+  +  H2, 

0(3P)  +  H2,  and  0(1D)  +  H2- 

•  3-D  translational  wavefunctions  have  been  obtained  for 
F  +  H2,  along  with  density  and  flux  maps.  The  3-D  flux 
maps  show  whirlpool  structure  similar  to  the  F  +  H2 
collinear  reaction  previously  studied.  This  will  be 
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presented  in  a  forthcoming  paper.  This  is  an  unusual 
method  of  interpreting  scattering  calculations  and  should 
lead  to  an  increased  understanding  of  molecular  reaction 
mechanisms. 

(14) 

e  An  integral  equation  method  developed  previously  ' 
has  been  tested  against  some  of  the  more  modern 
algorithms'*^ ,  and  in  many  instances  is  seen  to 
be  competitive.  This  algorithm  is  expected  to  be 
particularly  useful  in  applications  to  energy-dependent 
potentials  such  as  occur  in  reactive  scattering  problems. 

t  Several  scattering  codes  incorporating  electronic  coupling, 
using  different  integrators,  have  been  written  and  tested. 
Two  of  the  integrators  use  accurate  and  reliable  finite- 
difference  methods.  The  others  use  more  efficient 
potential-following  techniques.  The  finite  difference 
codes  can  be  used  to  test  the  accuracy  of  the  potential- 
following  codes  during  preliminary  studies  on  new  systems. 

•  Adiabatic  potential  energy  curves  and  nonadiabatic  first- 
derivative  couplings  for  the  X,  A,  and  C^E+  states  of  KH 
have  been  obtained  by  an  ab  initio  pseudopotential  method. 
The  important  splitting  between  the  X  and  A  curves  is  in 
good  agreement  with  experiment.  These  curves  and  couplings 
are  useful  for  dynamical  studies  on  this  system  (see 
Appendix  C). 
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•  The  ab  initio  potentials  were  used  to  calculate 
electronically  inelastic  transition  probabilities 

and  cross  sections  for  low-energy  K  +  H  collisions. 

2  2 

The  4  P  -*■  4  S  quenching  cross  section  varies  between 
2  x  10~4  a^  and  10  x  10~4  a^  between  .022  eV  and 
1.10  eV  relative  translational  energy.  This  study 
is  a  prelude  to  the  study  of  K  +  Ho. 

The  ultimate  goal  of  this  program,  namely,  the  treatment 
of  electronic  transitions  in  a  3-D  atom-diatomic  molecule  reaction, 
was  not  realized  due  to  time  constraints.  The  manner  in  which  the 
present  study  can  be  extended  to  this  process  is  discussed  in  the 
next  section. 

\T. _ RECOMMENDATIONS  FOR  FURTHER  WORK 

Progress  to  date  has  been  made  in  (1)  developing  efficient 
computational  tools  for  integrating  coupled  equations,  (2)  studying 
3-D  chemical  reactions  on  single,  adiabatic  potential  energy  surfaces, 
and  (3)  developing  a  formalism  for  including  electronic  transitions  in 
atomic  collisions.  The  following  recommendations  should  receive  serious 
consideration  to  fully  utilize  the  effort  expended  on  this  project. 

•  Perform  reactive  scattering  calculations  for  0(  P)  +  H^. 

Ab  initio  potential  surfaces  exist  for  this  system,  as 

do  classical  trajectory  results  and  experimental  results. 
This  will  provide  yet  another  reaction  for  which  3-D  quantum 
results  are  available.  This  was  scheduled  for  the  current 
project,  and  all  necessary  codes  are  in  hand. 
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•  Incorporate  electronic  degrees  of  freedom  into  existing 
scattering  codes.  This  was  started  during  this  project, 
and  the  3-D  reactive  code  has  been  partly  generalized. 

This  is  a  straightforward  extension  of  methods  currently 
in  hand,  at  least  for  singlet  diatomic  fragments. 

t  Develop  strategies  for  treating  arbitrary  electronic 

angular  momen*  im  in  a  total  angular  momentum  representation. 

•  Investigate  necessary  decoupling  approximations  suitable 
for  studying  vibronic  transitions. 

t  Perform  electronically  nonadiabatic  nonreactive  calculations 
on  Na  +  H2  or  K  +  H^.  This  will  test  the  machinery  necessary 
for  the  reactive  problem. 

•  Perform  electronically  nonadiabatic  reactive  calculations 
on  Na  +  H2  or  K  +  F^.  These  reactions  are  endothermic  by 
about  2.3  eV  and  2.7  eV,  respectively,  and  will  require  a 
large  number  of  channels.  Reliable  approximations  make 
these  problems  tractable  on  existing  machines. 

SCIENTIFIC  PERSONNEL  SUPPORTED  DURING  THIS  PROJECT 

Dr.  Bruce  C.  Garrett 

Dr.  Michael  J.  Redmon 


Dr.  Isaiah  Shavitt 
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Abstract 

Results  from  recent  three-dimensional  natural  coordinate  reactive  scattering  calculations  are 
presented.  Extensions  of  the  scattering  method  of  Wyatt  to  systems  with  nonlinear  intermediates 
are  discussed.  Rale  constants  for  the  reaction  H  +  H:  (t  =  1 1  at  300  K  are  presented  and  compared 
with  classical  trajectory  calculations  and  with  experiment  The  quantum  results  are  in  reasonable 
agreement  with  experiment,  but  the  classical  results  greatly  underestimate  the  reaction  rate  Totai 
cross  sections  and  relative  rate  constants  are  presented  for  the  F  +  H;  (r  =  0)  reaction  and  compared 
with  classical  results  and  experiment.  Total  cross  sections  for  the  H  +  O;  reaction  are  presented 
that  demonstrate  the  enhancement  of  reaction  caused  by  reagent  vibrational  energy 


1.  Introduction 

There  has  been  considerable  progress  in  the  development  of  quantum  me¬ 
chanical  methods  for  obtaining  state-to-state  cross  sections  and  rate  constants 
for  simple  chemical  reactions.  Beginning  with  the  early  work  on  the  H  +  H; 
reaction  [1-4),  the  computational  technology  has  continuously  developed  so 
that  converged  close-coupled  results  now  exist  for  this  system  [5-7],  Recently, 
a  calculation  has  been  reported  [7]  on  an  accurate  fit  [8]  to  the  definitive  Liu- 
Siegbahn  surface  [9],  It  should  now  be  possible  to  obtain  accurate  ab  initio 
dynamical  information  for  H  +  H;  and  its  isotopes  for  comparison  with  exper¬ 
iment. 

Progress  has  also  been  made  in  developing  methods  based  on  close-coupling 
techniques  for  treating  3-D  (three-dimensional)  systems  other  than  H  +  H; 
[  1 0- 1 2]  -  Extensions  to  reactions  involving  heavier  atoms  are  difficult  due  to 
asymmetries  in  the  reaction  coordinates  and  to  the  enormous  increase  in  the 
number  of  coupled  channels  necessary  for  convergence  of  the  computed 
probabilities.  The  large  number  of  channels  accessible  at  thermal  collision 
energies  requires  the  use  of  centrifugal  decoupling  approximations  for  total 
angular  momentum  J  >  0.  These  can  be  so-called  7r -conserving  approximations 
[13],  in  which  the  number  of  channels  used  in  expanding  the  wavefunction  is 
approximately  the  same  as  for  J  =  0.  or  centrifugal  sudden  approximations  [14], 
in  which  the  orientation  of  the  system  is  frozen  during  a  collision,  resulting  in 
an  enormous  reduction  in  the  number  of  coupled  equations.  Both  of  these  ap¬ 
proximations  are  adequate  for  total  reaction  cross  sections  for  H  +  H;.  and  the 
./.-conserving  approximation  reproduces  accuraie  close-coupling  results  for 
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individual  rotational  transitions  to  1%  1 10]  It  has  so  far  not  been  possible  to  test 
these  approximations  for  heavier  systems  by  comparison  with  close-coupling 
calculations.  It  is  felt  that  the  J. -conserving  approximation  should  be  reliable 
for  other  reactions  with  linear  intermediates. 

In  this  article  we  present  a  review  of  some  recent  results  of  3-D  reactive 
scattering  calculations  for  F  +  H;,  H  +  H;.  and  H  +  O;  These  calculations  used 
Wyatt's  formulation  of  the  3-D  quantum  scattering  problem  [15]  with  the 
natural  collision  coordinates  (see)  of  Marcus  [16).  The  calculations  employed 
the  code  REACTOR,  written  by  the  author  while  at  the  University  of  Texas 
at  Austin,  and  which  was  used  in  the  previous  calculations  reported  for  F  +  H; 
[10-12],  The  latest  version  incorporates  modifications  necessary  for  treating 
reactions  with  a  nonlinear  intermediate,  such  as  H  +  O:.  and  also  uses  the  /?- 
matrix  propagation  method  [17],  These  extensions  are  examined,  although 
computational  details  are  not  presented  here.  The  approximations  employed 
are  discussed  insofar  as  they  might  be  expected  to  affect  the  reported  results. 
A  recent  review  has  been  given  by  Wyatt  [18]. 

In  Sec.  2  we  discuss  some  improvements  in  the  ncc  approach  In  Sec  3  we 
present  results  for  several  systems,  including  rate  constants  for  the  F  +  H;  re¬ 
action.  and  make  comparisons  with  experiment  and  classical  trajectory  calcu¬ 
lations.  We  also  compare  our  calculations  of  reaction  rates  for  vibrationally  hot 
H  +  H;  with  some  new  experiments  and  with  classical  mechanics.  Finally,  we 
discuss  new  applications  of  the  method  to  the  H  +  0;  reaction. 

2.  Recent  Developments  in  NCC  Reactive  Scattering  Methodology 

A  More  Schizophrenia  in  Reaction  Coordinates 

The  unique  feature  of  Marcus'  natural  collision  coordinates  [16]  is  that  the 
translation,  vibration,  and  rotation  coordinates  (s,  p.  y)  all  vary  smoothly  from 
a  set  appropriate  for  describing  the  relative  motion  of  an  atom  C  colliding  with 
a  diatomic  molecule  AB  to  a  set  appropriate  for  describing  the  relative  motion 
of  atoms  A  or  B  with  the  molecules  BC  or  AC  This  is  accomplished  in  the  fol¬ 
lowing  way.  The  body-fixed  (.t.  A)  =  (0. 0)  plane  is  chosen  as  the  instantaneous 
plane  of  the  three  atoms,  with  coliinear  motion  defined  in  the  (r,  Z)  plane  De¬ 
viations  from  coliinear  motion  require  excursions  of  the  system  into  the  (v.  T) 
plane,  and  are  defined  by  the  quantity  m  >  0  with  magnitude  ( y:  +  }':) 1  :.  As 
discussed  by  Marcus  [16],  keeping  m  positive  avoids  one  source  of  double 
counting  of  configurations.  The  smooth  transformation  between  reactant  and 
product  coordinates  is  obtained  by  requiring  a  local  Cartesian  constraint 

y  sina(s)  +  5cos«U)  =  0  (I) 

at  each  point  along  the  reference  curve  defined  by  the  translational  coordinate 
s.  a(s)  is  an  arbitrary  switching  angle  [16]  that  varies  smoothly  between  zero 
for  reactant  configurations  and  £„  (skew  angle  of  the  mass-weighted  coordinates! 
for  product  configurations.  The  c  axis  points  initially  to  the  reagent  atom,  and 
switches  smoothly  so  that  it  points  to  the  product  atom  after  the  collision. 
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Rather  than  use  the  coordinate  m.  it  is  more  convenient  to  discuss  the  motion 
by  using  the  local  radial  and  bending  coordinates  (r.  7 ).  to  which  m  is  related 
by  the  expression 

m  =  r  siny;  (2) 

wnere  r  is  measured  (in  a  constant  s  plane)  from  the  r  axis  and  y  is  referenced 
to  the  collinear  plane.  This  is  useful  in  writing  the  Hamiltonian  for  the  system 
since  7  becomes  a  convenient  coordinate  for  representing  internal  rotational 
motion.  A  problem  that  arises,  which  was  not  discussed  by  Marcus,  is  that  there 
are  values  of  7  (near  ir/2)  for  which  the  coordinates  become  multiply  valued 
whenever  the  switching  angle  at(s)  is  not  exactly  equal  to  one  of  its  asymptotic 
values.  This  is  a  direct  consequence  of  the  constraint  expressed  by  Eq.  I 1 ).  This 
situation  is  illustrated  in  Figure  1.  where  the  hatched  area  for  7  >  y„  represents 
the  region  in  which  configurations  are  identical  to  some  for  7  <7*  Configu¬ 
rations  for  7  =  y„  correspond  to  isosceles  triangle  geometries,  and  when  s  = 
0  and  «  =  V;  {<*■  ym  -  -15°,  which  corresponds  to  equilateral  triangle  configu¬ 
rations  for  H3  and  H?+. 

This  characteristic  of  NCC  was  not  noted  in  early  applications  [6.  10.  11] 
because  the  bending  potential  was  parametrized  and  fit  to  small  deviations  from 
linearity  (small  7).  It  was  observed  by  this  author  during  attempts  to  accurately 
represent  the  F  +  H;  surface  in  NCC.  It  becomes  particularly  important  for 
systems  with  stable  nonlinear  intermediates  such  as  H?+  and  HO;.  In  fact,  for 
H+  +  Hi.  the  most  stable  configuration  follows  the  curve  7  *  y„  in  Figure  2. 
corresponding  to  C\,  configurations.  It  seems  necessary  to  replace  Eq.  <  2 ) 
with 

m  =  r  sin  by.  (?) 

with  the  scale  factor  b(s)  defined  so  that  67  =  ym  for  7  =  tr/2.  This  introduces 


Figure  1 .  A  constant  s  plane  showing  the  multivalued  region  for  7  >  y„  that  must  be  avoided 
in  doing  scattering  calculations.  When  the  switching  angle  equals  an  asymptotic  value.  1  „ 
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Figure  2.  A  comparison  of  ihe  variation  of  the  switching  angle  alj)  w  ith  reaction  coordinate, 
and  the  corresponding  value  of  ym  which  defines  isosceles  triangle  geometries  for  the 
three-atom  system 


additional  complications  into  an  already  formidable  kinetic  energy  operator  [15]. 
but  some  procedure  for  maintaining  7  <  y  is  necessary  for  accurate  compu¬ 
tations  in  NCC  It  should  be  pointed  out  that  this  schizoid  region  is  different  from 
the  one  associated  with  three-atom  dissociation  regions  of  NCC  [19).  which  are 
effectively  handled  for  low-energy  collisions  by  using  a  circular  arc  to  define 
the  reference  curve,  and  choosing  an  appropriate  turning  center. 


B.  Extensions  to  Systems  with  \onccllinear  Reaction  Paths 

NCC  theory  was  originally  formulated  with  applications  to  H  +  H;  in  mind 
Since  for  this  reaction  the  minimum  energy  path  is  collinear.  terms  in  the  kinetic 
energy  that  are  small  except  for  large  deviations  from  coliinearity  were  dropped, 
and  others  were  approximated  by  evaluating  them  on  the  reaction  path  [6.  I  5] 
This  near-linear  intermediate  approximation  has  been  used  in  all  see  calcu¬ 
lations  reported  so  far.  We  have  recently  added  the  additional  terms  that  con¬ 
tribute  for  J  =  0  and  now  neglect  only  those  that  are  zero  within  the  -/--con¬ 
serving  approximation.  Elkowitz  has  suggested  that  the  inertia  coefficient;,  can 
be  evaluated  on  the  noncollinear  reaction  path,  in  the  spirit  of  the  linear  inter¬ 
mediate  approximation  [20] .  He  has  shown  that  the  Hamiltonian  then  reduces 
exactly  to  the  one  used  previously  [  1 5]  when  the  reaction  path  is  collinear  We 
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have  chosen  instead  to  accurately  compute  the  matrix  elements  over  the  vibra¬ 
tional  motion,  which  can  be  large,  but  to  set  •>  equal  to  To-  its  local  value  on  the 
reaction  path,  to  approximate  some  of  the  integrals  over  the  bending  coordinate. 
This  leads  to  a  considerable  simplification  in  matrix  element  computation,  and 
seems  reasonable  since  we  are  looking  toward  developing  useful  decoupling 
procedures,  and  not  exact  close  coupling.  The  new  methods  of  evaluating  matrix 
elements  were  used  for  the  H  +  O;  calculations  discussed  in  this  article 


3.  Selected  Results  for  Three  Representative  Reactions 

In  this  section  we  present  some  recent  results  for  the  H  +  H;.  F  +  H;.  and 
H  +0:  reactions.  These  systems  are  useful  for  demonstrating  the  range  of  ap¬ 
plicability  of  the  method,  as  F  +  H2  is  highly  exoergic.  w  hile  H  +  O;  is  endoergic 
and  has  a  stable  nonlinear  intermediate  w  ith  a  2-eV  well.  The  goal  of  our  current 
effort  is  to  develop  techniques  for  determining  total  state-to-state  cross  sections 
and  rate  constants  for  atom-diatomic  molecule  reactions  involving  relatively 
light  atoms. 


-I.  Rate  Constants  for  the  H  +  Hz  Reaction 

In  their  work  on  this  reaction,  Schatz  and  Kupperman  [5]  computed  rate 
constants  for  the  Porter- Karpl us  surface  and  found  good  agreement  with 
classical  mechanics  at  600  K.  However,  at  300  K  there  were  significant  differ¬ 
ences  between  the  quantum  and  classical  results.  It  is  of  interest  to  compare 
quantum  and  classical  calculations  on  a  surface  with  the  correct  barrier  height, 
and  for  vibrationally  excited  reagents,  since  reactions  of  vibrationally  hot  hy¬ 
drogen  are  of  current  astrophvsical  interest.  We  have  chosen  the  Yates-  Lester 
surface  [21  ],  and  computed  distinguishable  atom  rate  constants  for  the  processes 
H  +  H2  (r  =  O.j  =  0)  —  H  +  H;  (c  =  0.  If),  and  H  +  H;  (;  =  )../  =  0)  -*•  H 
+  H2  (r  =  0, 2 j’).  The  results  are  summarized  in  Table  1.  We  find,  as  did  Schatz. 
that  the  classical  rates  at  300  K  are  significantly  lower  than  the  quantum  results. 
For  the  ground-state  reaction,  this  is  presumably  due  to  tunneling  For  vibra¬ 
tionally  excited  hydrogen,  tunneling  is  probably  less  important,  and  the  enhanced 
rate  is  the  result  of  interference  effects  among  the  various  reactive  and  non- 
reactive  pathways  that  suddenly  become  assessable  near  a  threshold. 

We  find  that  the  cooling  rate  for  H;  (r  =  1 )  is  slightly  larger  for  the  non¬ 
reactive  pathway  than  for  the  reactive  pathway,  as  predicted  from  J  =  0 
probabilities  [  1 8 j .  This  is  in  contradiction  to  the  assumptions  made  by  Heidner 
and  Kasper  in  analyzing  their  experiments  [22],  where  the  nonreactive  contri¬ 
bution  to  the  cooling  was  assumed  negligible. 

In  comparison  with  recent  hydrogen  maser  experiments  [23).  we  find  that 
our  overall  rates  are  generally  in  much  better  agreement  w  ith  experiment  than 
classical  results  on  the  Yates-Lester  surface  (at  300  K).  although  we  are  perhaps 
underestimating  the  contribution  due  to  the  resonant  exchange  process  H  +  H- 
U  -  1)-*H  +  H;(i  =  1 ).  W  e  are  currently  examining  the  possibility  that  the 
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Table  I  Rale  constants*  for  the  H  +  H;  reaction  at  300  k 


classical® 

£ 

Quantum 

experiment 

- 

1.3  X  lo8 

C1 .2  X  138 
d3.3  X  10s 

A 

5.1  *  ID10 

(3  X  1012) 

d3.1  X.  1C12 

A 

(2.4  x  to10; 

8.2  X  1CU 

- 

So 

- 

(1  X  1  o' 3 ) 

*1.8  X  1011 

‘!o 

3.6  X  1010 

2.6  *  10* 

- 

*  L  nits  are  cm3  sec**1  mole-1  Values  in  parenthesis  are  estimates 

”  A,  refers  to  reaction  from  H:  U  =  0.  j  «  0)  to  all  final  states,  A, o  refers  to  reaction  from  H: 
w  .  i  *  Oho  H;  <t  «  0.  Zj'i. 

J  Classical  results  on  the  Yates- Lester  surface  from  I  W  M  Smith.  Chem  Ph>v  Lett.  47.  21° 
<10*7) 

b  3-D  quantum  results  for  the  Yates- Lester  surface,  this  work 
c  \\  Schultz  and  D  J  LeRo\,  J.  Chem.  Phvs  42*  3869  U965) 
d  Reference  23. 
c  Reference  22. 

linear  intermediate  approximation  might  lead  to  underestimation  ot'  the 
probabilities  for  this  process.*  ti  should  be  noted  that  these  rates  are  computed 
from  distinguishable  atom  cross  sections  for  comparison  with  classical  me¬ 
chanics.  and  are  for  reagent  H;  in  its  lowest  rotational  state.  Our  conclusion  is 
that  the  use  of  classical  mechanics  for  this  system  is  justified  only  for  transla¬ 
tional  temperatures  well  above  300  K. 

B  Quantum  Effects  in  the  Three-Dimensional  F  +  Hz  Reaction 

This  reaction  was  the  first  one  studied  beyond  H  +  H;  by  a  full  3-D  quantum 
mechanical  method  (10.  1 I  j.  One  of  the  important  reasons  for  studying  this 
system  was  to  see  to  what  extent  the  very  significant  differences  between 
quantum  and  classical  collinear  calculations  (24)  might  be  modified  in  three 
dimensions.  The  original  3-D  quantum  calculation  was  restricted  to  total  angular 
momentum  J  =  0.  but  it  showed  that  the  Feshbacn  resonance  mechanism  that 
dominates  the  low-energy  collinear  reaction  probability  [25]  for  the  process  F 
+  H:  (r  =  0)  —  H  +  HF  <e'  =  2)  persists  in  three  dimensions.  It  was  later 
demonstrated  by  summing  over  J  to  obtain  a  total  cross  section  for  this  process 
that  the  quantum  cross  section  had  a  distinct  maximum  just  above  threshold, 
while  the  classical  result  continued  to  grow  [12], 

Our  current  quantum  scattering  codes  are  efficient  enough  to  allow  compu¬ 
tations  in  which  the  potential  parameters  of  the  surface  are  varied.  This  allows 
us  to  examine  the  sensitivity  of  the  dy  namics  to  various  features  of  a  potential 
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surface.  This  approach  can  be  used  to  aid  those  doing  ab  initio  calculations  of 
potential  energy  surfaces  in  selecting  those  regions  of  a  surface  that  require 
significant  effort.  This  has  been  done  successfully  for  the  collinear  F  +  H;  re¬ 
action  [25.  26],  Weare  presently  performing  3-D  calculations  for  F  +  H;  with 
a  variety  of  bending  potentials  and  present  in  Figure  3  cross  sections  for  a  fit  to 
the  famous  surface  5  of  Muckerman  [ 27 ] .  This  differs  from  the  potential  used 
previously  [10,  1 1  j  in  that  this  fit  gives  a  better  description  of  the  local  rotor 
eigenvalue  spectrum.  The  previous  potential  allow ed  for  less  hindered  rotational 
motion  near  the  saddle  point  than  the  one  used  here,  and  produced  reaction  cross 
sections  about  30°c  larger  near  the  threshold.  Results  for  all  of  the  potentials 
we  have  used  are  similar,  and  all  show  a  maximum  in  the  v  =  0  — ►  2  cross  section. 
As  Figure  3  shows,  this  produces  a  leveling  off  in  the  total  reaction  cross  section 
that  is  not  observed  in  the  classical  result  [12)  This  is  direct  evidence  that  a 
quantum  mechanical  resonance  mechanism  results  in  an  energy  dependence 
in  the  reaction  cross  section  that  is  not  reproduced  by  classical  mechanics,  which 
should  be  experimentally  verifiable. 

We  have  calculated  state-to-state  rate  constants,  and  present  them  in  Table 
II.  along  with  classical  and  experimental  results  where  possible.  We  find  that 
67 9b  of  the  available  energy  ends  up  as  product  vibration,  in  excellent  agreement 
with  experiment  and  classical  mechanics.  The  ratio  of  rate  constants  k,  k ;  is 
also  in  good  agreement,  but  the  ratio  k  - ,  k ;  is  not.  We  find  a  significant  amount 
of  flux  ends  up  in  r'  =  0  excitation  and  is  significantly  affected  by  variations  in 
the  bending  potential,  as  is  to  a  lesser  extent  v'  -  3.  Classical  mechanics  produces 
no  reaction  into  c'  =  0,  and  collinear  quantum  calculations  show  less  reaction 
into  this  state  than  3-D  calculations.  Modifications  to  the  bending  potential  can 
introduce  wells  in  the  c  -  0  and  r  =  1  correlation  diagrams  and  affect  the  amount 
of  reaction  into  those  states,  just  as  modifications  in  the  vibrational  levels 
themselves  can  drastically  affect  collinear  probabilities  [25]  This  sensitivity 
to  features  of  the  bending  potential  should  ultimately  lead  to  refinements  in  our 
knowledge  about  this  system. 


Figure  3  Cross  sections  Tor  the  reaction  F  +  H?  (r  =  0.  /  *  0>  -  H  +  HF  u\  £;'»  The  total 
reaction  cross  section  is  shown  and  compared  with  the  classical  result  of  ref  3>> 
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Table  II  Comparison  of  theoretical  and  experimental  results  for  the  3-D  F  +  H;  reaction 


classical 

.  t 
quantum 

experiment 

E,/kca1  mol*' 

1.9374 

l .  55 

I.6C 

log  A/ cm**  s‘  mol ' ' 

13. 2a 

13.6* 

14. 2C 

,f 1 
•  V 

.665- 

.c” 

,c£d 

k  ./k_  (3C3'>K'a,° 

v  max 

v'  »  0 

v‘  *  1 

.22 

.  24 

.31 

V1  »  2 

1.30 

1 . 02 

1.00 

v  *  3 

.26 

3 

.47 

0  ,/Q  {£  ,  «  2  ev 

v  max  v  rel 

v 1  *  0 

v  ’  *  1 

.23 

26 

.. 

v'  *  : 

1.30 

i  *0 

-- 

v 1  *  3 

.64 

.47 

-- 

‘  Classical  results  of  J  C.  Poianw  and  J  L  SchreiDer.  Faraday  Discuss.  Chem.  Soc.  62.  26" 
(1977);  T  and  ;  selected  from  a  300- K.  Boltrmann  distribution 
6  This  work.  Quantum  calculations  for  H;  If  =  0. )  -  0).  HF  t  Jr  .  -j'\:  thus  log  A  should  be 
underestimated. 

c  K.  H.  Homann  er  al  Ber  Bunsenges.  Phss  Chem  74.  fx?  1 19"0i. 
d  J  C.  Polanvi  and  K  B  Woodall.  J  Chem  Plus  5T.  1 5'4  1 19*2) 
e  Classical  cross  sections  of  ref  28 
r  Quantum  cross  sections,  this  work  and  ref  )  2. 

C.  A  Quantum  Mechanical  Study  of  the  H  +  0;  Reaction 

We  have  applied  the  nonlinear  intermediate  version  of  REACTOR  to  the  H 
+  O;  reaction,  using  a  LEPs  surface  of  Gauss  with  an  angle-dependent  Sato 
parameter  [29 1  This  reaction  is  of  practical  interest  because  it  is  an  important 
chain  propagation  step  in  many  combustion  systems  and  difficult  to  study  ex¬ 
perimentally.  Thus,  a  detailed  theoretical  investigation  of  the  kinetics  of  this 
system  is  warranted. 

From  a  theoretical  viewpoint,  this  reaction  presents  many  new  features.  It 
is  endothermic  by  about  V:  eV  and  has  a  metastable  intermediate  with  a  2-eV 
well.  The  integration  of  classical  trajectories  for  this  system  is  complicated  by 
the  occurrence  of  many  long-lived  complexes  which  often  lead  to  trajectories 
that  cannot  be  back-integrated  [ 29 ] .  The  quantum  calculations  are  made  dif¬ 
ficult  by  the  existence  of  so  many  open  channels. 

In  Table  III  we  present  total  cross  sections  for  formation  of  OH  (r  =0)  from 
various  vibrational  states  of  O;.  The  calculations  are  not  fully  converged  due 
to  core  limitations  on  the  CDC  7600  (we  were  limited  to  about  70  channels I. 
We  are  presently  using  a  VAX  1 1/780  where  core  size  is  not  a  limitation  and 
hope  to  converge  the  J- -conserving  calculations  for  this  system 

The  dependence  of  reaction  cross  section  on  reagent  vibrational  energy  follows 
the  trends  expected  for  an  endoergic  system  from  the  work  of  Polanyi  and  co¬ 
workers.  This  trend  was  also  noted  by  Gauss  for  this  system.  The  classical  cal- 
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Table  111.  Cross  sections  at  1-eV  total  energy  for  the  process  H  +  O;  l<  *  0.  /  «  Oi  -*  0  +  OH 
(«•'  *  0.  £/>. 


V 

Er.l 

3  „'A‘ 

3 

.503 

.15 

1 

.711 

.26 

2 

.522 

.30 

3 

.338 

.77 

£ 

.15? 

3 

culations  showed  no  tendency  of  O;  (r  =  0)  to  react,  yet  the  quantum  results 
indicate  thai  it  should  Gauss  essentially  found  no  reaction  below  t  =  4.  We  find 
the  cross  sections  are  small  for  r  <  4.  but  are  considerably  larger  than  the  clas¬ 
sical  result.  Gauss  estimates  that  his  t  =  4  cross  section  could  be  low  by  as  much 
as  a  factor  of  10.  due  to  the  inability  of  his  integrator  to  follow  many  complex 
tra  jectories.  We  are  presently  trying  to  estimate  the  error  in  our  result  due  to 
the  limited  basis  set  employed. 

We  find  that  nonreactive  collisions  tend  to  produce  vibrational!)  excited  O;. 
which  in  turn  can  react  rapidly  with  hydrogen.  A  detailed  study  of  this  system 
would  probably  produce  many  new  results  and  add  considerably  to  our  know  I- 
edgeof  elementary  combustion  processes. 
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Evidence  for  a  quantum  resonance  in  the  three-dimensional  F  +  Hj  (u  =  0,/  =  0)—  FH(u'  =  2,  all /')  +  H  reaction  is  pre¬ 
sented.  Relative  to  the  collinear  reaction,  this  resonance  is  much  broader  and  is  shifted  by  about  0.1  eV  to  higher  energies. 
This  resonance  has  not  been  predicted  in  previous  quasiclassical  trajectory  computations,  or  in  approximate  quantum  cal¬ 
culations. 


1 .  Introduction 

The  F  +  Hi  chemical  reaction  has  been  the  subject 
of  extensive  experimental  [1  ]  and  theoretical  study 
[2] ,  at  least  in  part  because  of  its  dominant  role  in 
powerful  F2/H2  chemical  lasers.  On  the  theoretical 
side,  the  reaction  has  been  attacked  from  the  view¬ 
points  of  classical,  semiclassical,  and  quantum  dynam¬ 
ics.  Prior  quantum  studies  consist  of  exact  collinear 
computations  on  several  different  potential  surfaces 
[2]  and  approximate  (Bom  [3]  and  distorted-wave 

[4] )  three-dimensional  treatments.  In  addition,  we 
have  reported  preliminary  quantum  results  on  the  3D 
reaction  which  are  based  upon  numerical  integration 
of  large  systems  of  quantum  close-coupled  equations 

[5] , 

In  this  study,  the  energy  dependence  of  reaction 
probabilities  and  cross  sections  for  the  reactions  F  + 
H2(u  =  0,/  =  0)  -»  FH(u',  all/')  +  H  are  reported  over 
the  total  energy  range  0.32  <£tot  <  0.50  eV 
=  £tot  -  0.27  eV).  Evidence  is  presented  for  a  broad 

*  This  research  was  supported  in  part  by  the  National  Science 
Foundation,  the  Robert  A.  Welch  Foundation,  and  Battelle 
Memorial  Institute. 


resonance  ir,  the  v  =  0  -»  v  -  2  cross  section,  with  a 
peak  in  the  cross  section  just  below  0.4  eV.  Resonance 
structure  of  this  type  has  not  been  predicted  in  either 
quasiclassical  trajectory  calculations  or  in  approximate 
quantal  calculations  [3,4],  However,  the  magnitude  of 
the  cross  section  and  the  energy  region  where  the  reso¬ 
nance  occurs  suggest  that  further  crossed  molecular 
beam  experiments  would  be  extremely  interesting  in 
testing  these  predictions.  A  brief  survey  of  the  scat¬ 
tering  methodology  is  presented  in  section  2.  New  re¬ 
sults  on  reaction  probability  surfaces  and  cross  sec¬ 
tions  are  then  presented  in  sections  3  and  4.  respec¬ 
tively. 

2.  Scattering  theory 

The  scattering  wavefunction  at  total  angular  momen¬ 
tum  J  is  expanded  in  products  of  adiabatic  hindered 
asymmetric  top  wavefunctions  [2,6],  12 ^  (80\y;  s), 
times  local  Morse  oscillator  functions,  HV(P'.  si 

*v0t0i0<.°*vrry m9 -112  S </0/0uj »(*) 
XS^*(0«xr.f)ffv(p:j),  (i) 
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where  {8,  o,  x}  are  Euler  angles  used  to  orient  the  top, 
and  j,  p,  y}  are  natural  translation-vibration-bending 
coordinates  [2,6],  Also  9  is  a  metric  coefficient  which 
is  used  to  scale  the  wavefunction  in  order  to  simplify 
the  structure  of  the  close-coupling  equations  for  the 
translational  wavefunctions,  ;  uy/(s).  In  the  com¬ 
putational  results  reported  here,  §(? rovibrational  chan¬ 
nels  were  employed,  with  the  distribution/ 12/1 2/ 1 2/ 
8/6/2/2,'2/2/2/,  where  the  total  number  of  even  j  plus 
odd  /  rotational  functions  in  each  of  the  ten  lowest  vi¬ 
brational  levels  is  indicated.  At  the  highest  energy  stud¬ 
ied  here,  six  of  these  vibrational  levels  are  asymptoti¬ 
cally  closed  in  products. 

For  all  computations  with  /  >  0.  the  J.  conserving 
approximation  [7]  was  employed  to  restrict  the  num¬ 
ber  of  orbital  angular  momentum  (/)  values  in  eq.  ( 1 ) 
to  a  single  “dominant"  term  for  each  value  of;  and  J. 
The  same  algorithm  for  selecting  l{j,  J)  was  success¬ 
fully  employed  in  earlier  H  +  H,  reactive  scattering 
calculations  [7],  In  that  case,  cross  sections  within 
about  17c  of  the  accurate  values  were  generated  with 
the  J,  conserving  approximation. 

The  close-coupled  equations  for  the  translational 
wavefunctions  were  numerically  integrated  with  the 
boundary  value  /{-matrix  propagation  method  [8]. 
Elements  of  the  5-matrix  were  then  directly  generated 
from  the  arrangement  channel  /{-matrices.  Reaction 
probabilities  for  F  +  H2(u  =  0./  =  0)  -*  FH(u',  all  /  ')  + 

H  reactive  collisions  are  defined  by 


=  /?o  /,£„  l5ooy— „yr<£)|2 


and  were  computed  for  total  energies  in  the  range  0.32 
eV  <  E  <  0.50  eV,  and  for  0  <  J  <  26  (in  ranges  of  J 
for  which  Pqu-  was  slowly  varying,  computations  were 
performed  at  every  other  J  value). 

In  these  scattering  calculations,  we  have  used  a  va¬ 
riety  of  potential  surfaces,  all  based  in  part  on  surface 
5  (M5)  of  Muckerman  [9]  *.  The  different  surfaces 
share  the  collinear  surface  of  MS,  but  they  differ  in 
the  range  and  degree  of  angular  anisotropy  of  the  bend¬ 
ing  potentials.  The  bending  potentials  all  have  the  same 
functional  form,  Kbend(7,  s)  =  5  ^(s)  (I  -  cos  2y), 
but  differ  in  the  position  along  the  reaction  coordinate 


(determined  by  s)  of  maximum  hindrance  to  rotation. 
The  surface  used  for  the  results  reported  here  is  simi¬ 
lar  to  the  M5  surface  for  small  deviations  from  colli- 
nearity  in  the  transition  state  region.  On  ihe  approach 
or  departure  from  the  transition  state,  the  bending  po¬ 
tential  is  similar  to  the  one  employed  in  recent  trajec¬ 
tory  studies  [12], 

3.  Reaction  probability  surfaces 

In  fig.  la  we  present  sections  through  the  reaction 
probability  surface  for  producing  HF(u'  =  2).  For  each 
value  of  the  total  angular  momentum  J.  the  variation 
of  probability  with  energy  is  qualitatively  similar  to 
that  obtained  from  a  collinear  calculation  on  the  M5 
surface  [10.13]  .  As  E  increases,  there  is  a  rapid  in¬ 

crease  and  then  decrease  in  the  probability,  followed 
by  a  fairly  constant  value  at  higher  energies.  One  strik¬ 
ing  effect  in  fig.  la  is  the  7-dependence  of  the  curves 
at  different  constant  energy  sections.  At  the  position 
of  the  /  =  0  resonance  maximum  and  below,  the  curves 
decrease  monotonically  with  J,  whereas  at  higher  ener¬ 
gies  the  curves  are  initially  fairly  constant  with  J.  but 
re-tune  onto  higher  resonance  values  at  a  value  of  J 
that  increases  with  energy.  The  locus  of  these  maxima 
in  the  E  -  J  plane  follows  the  “resonance  ridge”  that 
begins  for  J  =  0  at  the  maximum  in  the  probability 
curve,  and  progressively  moves  to  higher  J  with  in¬ 
creasing  energy.  The  maximum  probability  on  the  ridge 
gradually  decreases  as  E  and  J  increase  The  ^  curves 
for  the  3D  reaction  are  broader  than  the  collinear  re¬ 
sult  [10,13,14]  due  to  the  participation  of  many 
(twelve  in  these  calculations)  rotor  states  in  the  reso¬ 
nance  mechanism,  and  to  changes  in  the  shape  of  the 
vibrational-rotation  adiabatic  energy  correlation  curves 
with  increasing  values  of  J. 

In  our  earlier  collinear  calculations  on  this  reaction 
[  1 3  ] ,  we  have  identified  the  resonance  mechanism  in 
the  0  -*  2  process.  It  arises  from  interna)  excitation 
and  then  de-excitation  of  the  FHH  intermediate  in  the 
beginning  of  the  FH  +  H  exit  valley  into  primarily  the 
u  =  3  and  4  states  (which  may  be  approximately  iden¬ 
tified  as  asymmetric  stretch  states)  of  the  vibrational 
energy  correlation  curves.  Qualitative  similarity  be¬ 
tween  the  collinear  and  3D  results  (at  each  J)  suggests 

*• 


See  ref.  [10]  for  M5  parameters.  For  early  classical  results 
on  predecessors  of  the  MS  surface  see  ref.  [11]. 


For  extensive  bibliographies  see  refs.  [1 3.14). 
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Fig.  1.  Reaction  probability  turf  ace  for  F  +  Hj(u  =  0,  /  =  01 
—  H  +  HF  (u  -  2./'  3  all),  (b)  Reaction  probability  surface 
for  F  +  Hj(u  *  0,/  =  0)  —  H  +  HF(i>  =  3,/'  =  all).  Transla¬ 
tional  energy  =  total  energy  -  0.27  eV. 

that  the  same  mechanism  is  operating  in  the  3D  colli¬ 
sions. 

It  is  known  from  other  work  that  variations  in  the 
collinear  surface  in  the  saddle-point  and  downhill  re¬ 
gions  [13]  and  variations  of  the  overall  exothermicity 
and  product  vibrational  spacing  [14]  drastically  affect 
computed  reaction  probabilities,  by  shifting  or  elimi¬ 


nating  resonances  or  by  creating  new  ones.  The  effects 
we  observe  by  varying  the  bending  potential  are  not 
nearly  as  severe  as  those  resulting  from  changes  in  the 
collinear  surface,  and  mainly  serve  to  alter  the  magni¬ 
tudes  of  the  resonance  probabilities  and  the  values  of 
J  that  contribute  most  to  the  integral  cross  sections. 

4.  Reaction  cross  sections 

In  fig.  2  we  present  the  energy  dependence  of  cross 
sections  for  reaction  of  ground  state  H-,  to  form  HF. 
Also  shown  is  a  classical  total  cross  section  on  a  slight¬ 
ly  different  potential  surface  [12].  There  is  good  agree¬ 
ment  below  the  resonance  maximum,  but  above  that 
energy  the  quantal  and  classical  results  differ  substan¬ 
tially.  The  classical  result  continues  to  rise,  while  the 
quantum  total  cross  section  levels  off,  due  to  the  reso¬ 
nance  in  the  0  -+  2  cross  section.  In  earlier  collinear 
studies  [10],  the  classical  total  reaction  probability 
was  also  found  to  exceed  (by  about  a  factor  of  two) 
the  quantum  result,  between  the  classical  threshold  at 
0.29  eV  and  the  energy  region  where  the  0  -*  3  reac¬ 
tion  probability  begins  to  grow  (0.4  eV).  However,  the 
0  -*■  2  resonance  width  in  fig.  2  is  much  broader  than 
in  the  collinear  calculations  (where  the  width  is  about 


Fig.  2.  Cross  sections  for  forming  vibrationallv  excited  HF. 
for  all  open  HF  vibrational  manifolds.  The  classical  result  is 
from  ref.  [  1 2  J . 
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0.01  eVl.  The  broadening  mechanism  in  the  3D  case 
arises  both  from  the  participation  of  many  rotor  states, 
and  the  detailed  angular  momentum  dependence  of 
the  "resonance  ridge”  in  fig.  la.  In  our  studies  on  other 
potential  surfaces  with  different  bonding  potentials, 
the  0  —  2  cross  section  above  0.40  eV  never  showed 
smooth  growth  that  could  be  extrapolated  from  the 
lower  energy  region,  as  in  the  classical  3D  calculations 
[4-12],  The  quantum  0  —  2  cross  sections  for  differ¬ 
ent  bending  potentials  ranged  from  about  3.0  to  5.5 
Up  at  0.50  eV  and  thus  showed  either  a  definite  level¬ 
ing-off  or  a  slight  decline,  relative  to  the  values  near 
0.40  eV.  The  peak  0  —  2  cross  section  near  0.40  eV 
ranged  from  4  to  6  Jq.  except  for  one  model  bending 
potential  which  maximized  the  angular  anisotropy  near 
the  beginning  of  the  exit  valley  .  where  the  peak  cross 
section  was  1 1  ajj. 

Further  analysis  of  these  cross  sections  indicates 
that  brrc  of  the  available  energy  appears  as  product 
vibration  at  the  0  —  2  resonance  maximum  (0.4  eV). 
in  agreement  with  both  experiment  [1,15]  and  trajec¬ 
tory  calculations  [9-1 1  ]  in  the  low  energy  region. 

5.  Conclusions 

These  quantum  scattering  calculations  of  reaction 
probabilities  and  cross  sections  for  the  three-dimen¬ 
sional  low  energy  electronically  adiabatic  F  +  Hi  reac¬ 
tion  have  demonstrated  resonance  structure  for  the 
v  =  0  —  v  =  2  process  which  is  much  broader  and  at 
0.1  eV  higher  energy  than  the  collinear  resonance,  and 
which  is  not  predicted  by  quasiclassica!  trajectory  cal¬ 
culations.  Contribution  of  this  resonance  to  the  shape 
of  the  total  cross  section  has  not  been  shown  experi¬ 
mentally.  The  extension  of  previous  crossed  molecular 


beam  studies  (on  the  F  +  D-,  reactiun  |16|  im  the  reso¬ 
nance  energy  region  between  I'-Uans  =  0.05  eV  to/.  =0.23 
eV  could  test  these  predictions.  Resonance  structure  m 
this  reaction  should  be  more  amenable  to  experimen¬ 
tal  study  than  the  vibrational  resonance  predicted  for 
the  H  +  H-,  reaction  [17], 
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ABSTRACT 


The  accurate  ab  initio  MBPT  quartic  force  field  of  Bartlett, 

Shavitt  and  Purvis  has  been  fit  to  an  analytical  function  using  a  method 
developed  by  Sorbie  and  Murrell  (SM).  An  analysis  of  this  surface  indicates 
that  it  describes  most  properties  of  the  FLO  molecule  very  accurately, 
including  an  exact  fit  to  the  MBPT  force  field,  and  very  close  to  the  correct 
energy  difference  between  linear  and  equilibrium  H20.  The  surface  also  re¬ 
produces  the  correct  diatomic  potentials  in  all  dissociate  /e  regions,  but 
some  aspects  of  it  in  the  "near  asymptotic"  0(1D)  +  FI,  region  are  not  quan¬ 
titatively  described.  For  example,  the  potential  seems  to  be  too  attractive 
at  long  range  for  0  +  H„  encounters,  although  it  does  have  the  correct  minimum 
energy  path  geometry  and  correctly  exhibits  no  barrier  to  0  atom  insertion. 
Comparisons  of  this  surface  with  one  previously  developed  by  SM  indicates 
generally  good  agreement  between  the  two,  especially  after  some  of  the  SM 
parameters  were  corrected,  using  a  numerical  differentiation  algorithm  to 
evaluate  them.  A  surface  developed  by  Schinke  and  Lester  (SL)  is  more  realistic 
than  ours  in  the  0(:D)  +  FI,  regions,  but  less  quantitative  in  its  description 
of  the  H,0  molecule.  Overall,  the  present  fit  appears  to  be  both  realistic 
and  quantitative  for  energy  displacements  up  to  3-4eVfrom  H.O  equilibrium, 
and  should  therefore  be  useful  for  spectroscopic  and  collision  dynamics  studies 
involving  FLO. 


I.  INTRODUCTION 


In  a  recent  paper,  Bartlett,  Shavitt  and  Purvis'1  presented  the 
results  of  an  accurate  ab  initio  calculation  of  the  ground  state  quartic 
force  field  of  H.,0.  This  calculation  used  a  many  body  perturbation  theory 
(MBPT)  method  including  up  to  quadruple  excitations  and  a  large  basis  set 
(39-STO)  wave  function  to  evaluate  all  of  the  31  quadratic,  cubic  and 
quartic  force  constants  in  the  generalized  valence  force  field.  Not  all  o* 
these  force  constants  have  been  determined  experimentally,  but  where  accurate 
values  are  known,  the  MBPT  values  are  in  good  agreement  with  them.  Indeed 
the  MBPT  force  field  may  be  better  than  experiment,  but  in  order  to  use  this 
force  field  for  spectroscopic  or  scattering  calculations,  it  must  be  extended 
to  map  out  regions  of  nuclear  configuration  space  away  from  the  H,0  equilib¬ 
rium  geometry. 

2 

In  this  paper,  we  use  the  method  of  Sorbie  and  Murrell  to  fit  the 
MEPT  surface  to  an  analytical  function.  This  function  identically  reproduces 
the  MBPT  H20  quartic  force  field,  and  correctly  describes  the  0(*D)  +  H, 
and  Q'n(2p]  )  +  H  dissociative  channels  at  infinite  separation.  In  between 
these  limits,  a  smooth  interpolation  is  Drovided.  We  also  describe  a  simple 
numerical  algorithm  for  generating  the  parameters  used  in  the  Scrbie-Murrel 1 
(SM)  fitting  method,  and  we  correct  their  fit  to  the  spectroscopically 

3 

derived  force  field  of  Hoy,  Mills  and  Strey  (HMS). 

Besides  the  SM  surface,  other  complete  surfaces  for  H:0(1A1)  have 
been  generated  by  Tully^3,  by  Whitlock,  Muckerman  and  Fischer^  (WMF)  and  by 
Schinke  and  Lester^  (SL).  Tully  and  WMF  used  the  valence  bond  diatomics  in 
molecules  (DIM)  method  to  construct  their  potential  surfaces.  The  two 


m 
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surfaces  are  not  identical  because  of  differences  in  integral  evaluation, 
but  neither  surface  describes  the  H20  molecule  with  anything  close  to 
spectroscopic  accuracy.  SL  used  a  Sorbie-Murrell-like  function  to  fit  the 
ab  initio  surface  of  Howard,  McLean  and  Lester^  (HML),  but  unlike  SM  (see 
Section  II),  they  obtained  the  coefficients  in  their  fit  by  a  least  squares 
analysis.  A  fairly  widely  scattered  set  of  ab  initio  ooints  was  used  for  the 
rit,  with  the  result  that  their  surface  describes  the  0  +  H2  and  OH  +  H 
regions  more  accurately  than  SM,  but  their  H20  force  field  is  much  less 
accurate  (though  still  better  than  the  DIM  ones). 

Because  calculated  rather  than  experimental  dissociation  energies 
were  used  for  all  or  parts  of  the  potential  surfaces  or  Tully,  WMF  and  SL, 
certain  energetic  aspects  of  these  surfaces  are  in  error.  For  example,  the 
energy  associated  with  dissociation  of  K20  to  0(3P)  +  2 H ( 2 S )  is  10.08  eV 
experimentally,  but  only  9.37  eV  on  Tully's  surface,  9.16  eV  on  WMF's  and 
8.95  eV  on  SL.  Likewise,  the  0H(zn  )  dissociation  energy  is  4.63  eV  experi¬ 
mentally,  but  4.54  eV  on  Tully,  4.58  eV  on  WMF  and  4.04  eV  on  SL.  Both  Tullv 
and  WMF  do  correctly  describe  the  energy  associated  with  0(]D)  +  H,  ->-0(3P)  + 

2H ( 2 S )  (2.79  eV)  but  SL's  value  is  1.87  eV.  In  this  paper,  we  use  the  experi¬ 
mental  dissociation  energies  so  as  to  insure  the  proDe>'  energetics  in  all 
arrangement  channels. 

Because  the  ab  initio  MBPT  points  are  clustered  close  to  the  H20 
equilibrium,  the  present  fit  describes  the  H20  molecule  properties  much  more 
accurately  than  the  0  +  H2  and  OH  +  H  regions.  This  is  like  the  SM  surface, 
but  in  the  present  case,  certain  MBPT  derived  quartic  force  constants  not 
available  to  SM  have  been  included  in  the  fitting  process.  We  will  examine 
here  the  influence  of  these  additional  constants  on  the  long  range  nature  of 
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the  surface.  In  addition,  a  general  comparison  of  the  H?0  molecule  pro- 
oerties  of  the  SM,  M8PT  and  SL  surfaces  will  be  made. 

II.  EVALUATION  OF  FITTED  SURFACE 

2 

The  functional  form  of  the  fitting  surface  for  H„0  is 

vM(Ri>R2>R3)  »  VQH( Ri)  +  Vqh(R2)  +  Vkh(R3) 

+  V3(R1,R2,R3)  ,  (1 ) 

where  R3  is  one  of  the  OH  distances,  R2  the  other  and  R3  the  HH  distance. 

and  VHH  are  the  OH  and  HH  ground  state  potential  curves,  and  are  taker 
from  Ref.  7,  while  V3  accounts  for  all  3-body  terms  in  Q.  Not  ■  that 
Rj  <  1.6707A,  H.,0  dissociates  into  0(:D)  +  H,  (1J+)  but  for  R3  >  1.6707A,  tre 

X 

ground  state  dissociation  channel  is  0(3P)  +  "2  <3U-  VHH  incorporates  tv's 

i„  +  3„  + 

behavior  by  switching  from  the  )_  to  the  potential  curves  at  that  P-. 

value.  This  leads  to  a  cusp  in  both  V3  and  VH  g  asymptotically  and  neces- 

2 

sitates  the  use  of  a  discontinuous  form  for  V3  in  order  to  make  g  contin¬ 
uous  near  the  H„0  equilibrium  geometry. 

2 

The  functional  form  of  V3  is  given  by  : 

V3  =  A(l-tanhv1S1/2)(l-tanhY2s2/2)(l-tanhY3S3/2)P(S1,S2,S3)  (2) 

where  Si  =  R.  -  R°  (i  =  1,2,3)  and  R?  is  the  ith  H,0  equilibrium  internuclear 
distance.  A,  Ylt  y2  and  y3  are  parameters  and  P  is  a  polynomial  in  $lt  S; 
and  S3.  This  polynomial  will  have  a  different  representation  inside  the 
above  mentioned  cusp  than  outside,  and  we  denote  these  two  polynomials  as 
Pin  and  Pout,  respectively. 

As  discussed  in  Ref.  2,  the  coefficients  in  P1n  are  conveniently 
evaluated  when  an  analytical  representation  of  the  H.,0  force  field  is  known 


(near  equilibrium)  by  relating  the  derivatives  of  this  force  field  (evaluated 
at  equilibrium)  to  those  of  Pln.  Although  explicit  analytical  expressions 
for  these  polynomial  coefficients  were  given  in  Ref.  2,  their  evaluation  is 
tedious  for  all  but  the  simplest  types  of  force  fields.  A  much  easier 
evaluation  of  these  coefficients  can  be  accomplished  by  numerical  differen¬ 
tiation  of 


d  = 


/H-0'v0H^R’^"V0H^R2^“VHH^R3) 


(3) 


(l-tanhv:S1/2) (l-tanhv2S2/2)(l-tanhv3S3/2) 
ir.'hen  derivatives  of  (3)  are  evaluated  at  the  equilibrium  position,  the 
resulting  values  are  sirndy  proportional  to  the  polynomial  coefficients  in 
r',n.  P0Ut  can  be  similarly  evaluated  by  requiring  that  the  derivatives  of 
the  ootential  at  S3  =  S3  (cusp)  (with  S3  =  S2  =  3)  be  continuous.  By  using 
the  previously  obtained  ?ln  to  evaluate  V3  for  S3  slightly  inside  the  cuso 
in  Eq.  (3)  and  for  S3  slightly  outside,  numerical  differentiation  of  Eq. 
(3)  directly  yields  the  coefficients  of  the  Dolynomial  outside,  as  expanded 
ab.ut  the  cusp  position,  Pout  can  then  be  reexDanded  about  the  equilibrium 
Dosition,  if  desired,  by  evaluating  its  numerical  derivatives  at  that 
position.  Note  that  the  same  numerical  differentiation  program  is  used  three 
times  in  this  evaluation,  once  in  determining  Pin  and  twice  for  Pou  .  By 
orogramming  this  algorithm  in  double  precision  (64  bit  words)  and  using  a 
judicious  choice  of  finite  difference  increment  (1  x  10  bohr  for  the  first 
and  second  derivatives,  5  x  10~^  for  3rd  and  4th),  even  the  simplest  differ- 

g 

entiation  formulas  enable  the  determination  of  4th  derivatives  to  3-4 


significant  figures  (with  much  higher  significance  for  the  lower  order 
derivatives).  Moreover,  this  algorithm  is  independent  of  the  functional 
form  used  to  represent  the  a b  initio  potential  near  the  equilibrium  position. 
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The-  matching  procedure  at  S3  =  S3  (cusp)  does  not  guarantee  con¬ 
tinuity  of  the  potential  across  the  cusp  for  Si  t  S- .  This  deficiency  was 
corrected  by  switching  between  polynomials  at  the  cusd  using  the  expression 


where 


n  .innin  .  2  .  „outnout  : 

P  =  A  P  sin  u  +  A  P  cos 


a)  =  4  rl-tanhi'sCS2-S2 (cusp)]} 


(4) 

(5) 


The  value  of  Darameter  r  determines  the  range  of  mixing  of  the  two  poly¬ 
nomials  on  either  side  of  the  cusp. 

The  resulting  polynomial  coefficients  using  the  SDQ-M8PT(4) 
potential  force  field  of  Ref.  1  are  listed  in  Table  I  (labelled  MBPT).  Also 
given  are  the  other  parameters  in  V3  mentioned  previously  (rl5  v2  and 
are  taken  from  Ref.  2),  and  the  analogous  coefficients  and  parameters  in  a 
fit  (using  our  method)  to  the  soectroscopically  derived  force  field  of  Hoy, 
Mills  and  Strey3  (labelled  HMS).  This  latter  force  field  was  also  fit  by 
Sorbie  and  Murrell,  and  in  Table  I,  we  list  their-  parameters  (labelled  SM). 
Since  the  same  method  of  fitting  was  used  in  each  aco" 1 cation ,  the  SM  and 
HMS  Darameters  should  be  identical.  Table  I  indicates  tnat  most  of  the 
parameters  agree  to  3-4  significant  figures.  Tw~  sets  of  coefficients  are 
very  different  however.  The  SM  coefficient  multiplying  Si  is  of  similar 

magnitude  but  opposite  in  sign  to  ours  for  both  p1n  and  pout,  while  the 

0 

coefficient  multiplying  S 1 S2  differs  by  a  factor  of  about  3.  The  origin  of 
these  differences  is  not  known,  but  in  the  next  section,  we  shall  see  that 
they  cause  SM’s  potential  to  have  slightly  different  quartic  force  field 
parameters  than  HMS,  even  though  SM  is  supposed  to  be  a  fit  to  HMS. 

Comparison  of  the  HMS  and  MBPT  parameters  in  Table  I  indicates 
good  agreement  for  the  lower  order  coefficients,  but  only  qualitative  agreement 
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for  the  higher  order  ones.  The  influence  of  these  higher  order  coefficients 
on  features  of  the  potential  surface  will  be  considered  in  the  next  two 
sections. 


III.  PROPERTIES  OF  FITTED  SURFACES 
A.  H.O  Force  Field  Parameters 


The  first  requirement  of  the  MBPT  fitted  surxace  is  that  it  should 

be  able  to  reproduce  the  MBPT  quartic  force  field  exactly.  A  convenient 

representation  of  this  surface  involves  an  expansion  in  scaled  normal  mode 
9 

coordinates.  If  Q$  is  the  normal  coordinate  for  mode  s,  then  the  dimen¬ 
sionless  coordinate  q$  is  defined  by 


.  _  o  t Cus \ 1 / 2  n 

qs  “  2’  °s 


and  the  potential  VH  q  is  expanded  as 


'H„0 


/he  =  +  u2q22  +  w3q32) 


(6) 


+  kmq!3  +  k122q1o22  +  k133q1a;,:-  ~  k2nq2q:2 

+  k222q2  -  +  k2  33q2q32  +  knnq^  +  k1122q12a:: 

+  k3 1 33di2q  32  +  k2222q2  4  +  k22  33q22q32  +  k3333q3“  (7) 

Values  of  the  coefficients  in  this  expansion  for  the  fitted  MBPT  surface  are 
given  in  Table  II.  The  MBPT  fitted  surface  is  labelled  MBPT-F,  while  the 
quartic  representation  used  to  generate  the  fit  is  labelled  MBPT-Q.  Like¬ 
wise,  our  fit  to  the  HMS-Q  surface  is  labelled  HMS-F,  while  Sorbie  and 
Murrell's  fit  is  labelled  SM-F.  Also  included  in  this  Table  are  the  analogous 

5 

coefficients  obtained  from  SL's  surface. 
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We  note  first  in  Table  II  that  all  the  quartic  parameters  obtained 
from  our  fits  to  HMS-Q  and  MBPT-Q  agree  with  those  of  the  surface  being  fit 
to  within  0.3  cm”"'  or  better.  The  agreement  should,  of  course,  be  exact,  and 
for  most  parameters,  it  is.  A  few  parameters  are  slightly  different,  probably 
because  of  round-off  errors  in  the  4th  derivative  evaluation  used  to  generate 
the  fitted  surfaces. 

SM-F's  parameters  differ  from  HMS-0  by  as  much  as  7.7  cm”"',  probably 
because  of  the  differences  in  coefficients  noted  in  Table  I.  SL's  parameters 
in  Table  II  are  very  different  from  either  MBPT-Q  or  HMS-Q  (over  330  cm”"'  for 
u. )  although  most  of  the  parameters  have  at  least  Qualitatively  reasonable 
values.  The  large  differences  between  the  SL  results  and  those  of  the  other 
surfaces  are  at  least  partially  due  to  the  lower  accuracy  of  the  ab  initio 
results  being  fit,6  and  partially  to  the  least  squares  procedure  used  by  SL  to 
fit  HML's  points  (SL's  fitting  method  places  less  emphasis  on  the  H20  equilib¬ 
rium  geometry).  It  should  be  noted  that  the  equilibrium  H20  geometry  o*  these 
surfaces  shows  variations  similar  to  those  between  the  force  constants.  This 
is  indicated  in  Table  III,  where  the  OH  equilibrium  distances  and  H20  ancle 
are  listed.  Also  included  in  the  table  are  the  values  of  q  at  equilibrium 
for  a  number  of  H20  potential  surfaces.  As  noted  in  the  introduction,  the 
values  for  the  SL,  WMF  and  Tully  surfaces  differ  by  roughly  1  eV  from  the 
experimental  value. 

B.  Other  Properties  of  MBPT-F  and  Other  Surfaces 

In  Fig.  1  is  plotted  q  as  a  function  of  the  H20  bending  angle 
3  (with  Si  =  S2  =  0)  for  the  MBPT-Q,  MBPT-F  and  SL  surfaces.  The  correspond¬ 
ing  curves  for  HMS-Q  and  HMS-F  are  quite  similar  to  MBPT-Q  and  MBPT-F, 


8 


resoectively ,  in  the  figure.  Notice  how  the  MBP7-F  curve  is  considerably 
more  repulsive  than  MBPT-Q  for  small  3  and  less  for  large  £.  The  SL  curve 
is  qualitatively  similar  to  MBPT-F  for  large  9  but  not  for  small. 

Of  particular  interest  in  the  analysis  of  the  bending  potential 
is  the  characterization  of  linear  H20.  None  of  the  quartic  force  fields  are 
even  qualitatively  correct  for  that  geometry  since  the  correct  potential  has 
a  saddle  point  there  while  the  quartic  representations  have  a  nonzero  value 
of  iVH  g/39  at  3  =  ’7.  All  of  the  fitted  surfaces  do  have  saddle  points,  and 
the  oroperties  of  these  are  summarized  in  Table  IV.  Of  Darticuiar  importance 
in  this  table  is  the  energy  difference  PV  between  the  saddle  point  energy  and 
the  water  equilibrium  energy.  Experimental  estimates  of  this  barrier10 
indicate  a  value  1.37  eV,  which  is  in  best  agreement  with  the  MBPT-F  value  of 
1.296  eV.  The  HMS-F  value  (1.168  eV),  which  is  very  close  to  SM-F  (1.175  eV), 
is  too  low  while  SL 1 s  value  (1.708  eV)  is  much  too  high.  It  is  also  inter¬ 
esting  to  note  that  although  the  HMS-F  and  MBPT-F  R]  values  in  Table  IV  are 
only  slightly  smaller  ( <0.01  A)  than  at  equilibrium,  SL's  value  is  0.04A 
smaller.  SL's  saddle  point  frequencies  are  also  aporeciably  different  (300- 
300  cm-1)  from  MBPT-F  or  HMS-F. 

In  Figs.  2  and  3  are  plotted  contour  diagrams  of  the  MBPT-F  surface 
as  a  function  of  the  two  OH  distances  Rj  and  R2  for  the  equilibrium  and  linear 
H20  configurations.  A  cut  through  Fig.  2  corresoonding  to  symmetric  bond 
stretching  displacements  is  given  in  Fig.  4.  Included  in  this  figure  is  a 
comparison  of  the  MBPT-Q,  MBPT-F  and  SL  surfaces.  There  is  good  agreement 
between  HMS-F  and  MBPT-F  for  symmetric  bond  stretch  displacements  exceot  at 
very  large  and  small  values  of  Ri ,  so  we  have  not  plotted  the  former  curve  in 
the  figure.  Fig.  4  does  indicate  that  the  MBPT-Q  potential  is  more  repulsive 
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than  MBPT-F  at  both  large  and  small  Ri«  The  difference  between  MBPT-Q  and 
MBPT-F  remains  smaller  than  0.1  eV  for  displacements  over  2  eV  away  from 
equilibrium.  This  contrasts  with  the  behavior  in  Fig.  1,  which  indicates 
0.1  eV  differences  only  0.7  eV  above  equilibrium.  Since  the  analogous  com- 
Darison  for  the  asymmetric  bond  stretching  potential  indicates  excellent 
agreement  between  MBPT-Q  and  MBPT-F  for  several  eV  displacements,  we  conclude 
that  the  most  serious  errors  in  the  quartic  representation  of  the  H,Q  DOtentiai 
arise  for  bending  displacements.  These  errors  can  be  important  even  for  the 
low  lying  vibrational  states  of  H,0  since  the  vibrational  zero  Doint  energy 
(0.59  ev/ )  is  comparable  to  the  energy  displacement  needed  to  make  the  differ¬ 
ences  between  MBPT-F  and  MBPT-Q  large. 

Comparison  of  the  SL  and  MBPT-F  curves  in  Fig.  4  indicates  reasonable 
correspondence  for  large  Ri  corresponding  to  3-body  dissociation.  This 
indicates  that  at  least  for  coordinate  displacements  of  the  type  indicated, 
the  fitted  MBPT  surface  agrees  with  one  known  to  be  more  accurate  asymDtoti cally. 

C.  Surface  Properties  in  the  0(2D)  +  H,  Region 

In  Fig.  5  are  plotted  contours  of  the  M3PT-F  surface  corresoonding 
to  C2v  symmetry  collisions  of  0(1D)  with  Ho.  The  distance  X  in  the  figure 
is  the  oxygen  atom  to  center  of  mass  of  H2  distance.  Previous  ab  initio 
studies  have  indicated^  ^  "*  that  the  minimum  energy  path  for  0  (:D)  +  OH  +  H 

follows  this  and  similar  geometries,  and  the  present  MBPT-F  surface  concurs  with 

\ 

this.  Two  cuts  through  the  contours  in  Fig.  5  are  Dlotted  in  Figs.  6  and  7. 

Fig.  6  plots  the  MBPT-F,  HMS-F  and  SL  values  of  VH  g  versus  X  for  the 
equilibrium  H2  value  of  R3,  while  Fig.  7  presents  the  analogous  plot  for  the 
eouilibrium  H,0  value  of  R3.  No  barriers  are  evident  in  the  potential  curves 
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for  any  of  the  three  surfaces.  This  agrees  with  the  results  of  Ref.  5,  but 

disagrees  with  the  barrier  found  by  Gangi  and  Bader"'"'.  Since  the  experi- 

1 2 

mental  activation  energy  for  0( *D)  +  H2  is  apparently  zero  ,  we  presume  that 
the  zero  barrier  result  is  correct. 

Fig.  6  indicates  that  for  R3  equal  to  the  H,  equilibrium  value,  tne 
MBPT-F  and  HMS-F  potentials  exhibit  strong  attractive  behavior  at  much  larger 
X  values  than  SL.  Since  the  minimum  energy  paths  for  large  X  follow  the  curves 
in  this  figure,  we  conclude  that  the  HMS-F  and  MBPT-F  minimum  energy  paths 
have  much  more  attractive  profiles  for  large  X  than  SL.  Close  to  H20 
equilibrium,  the  minimum  energy  paths  switch  to  the  curves  depicted  in  Fig.  7, 
where  the  MBPT-F,  HMS-F  and  SL  curves  are  quite  similar.  Thus  it  is  only  at 
larger  distances  where  there  is  much  discrepancy  between  the  surfaces,  and  we 
find  that  even  in  this  limit,  the  MBPT-F  and  HMS-F  are  still  very  close.  This 
suggests  that  the  quartic  terms  in  the  force  fields  (which  differ  substantially 
between  MBPT-F  and  HMS-F)  do  not  play  an  important  role  at  the  larger  distances 
considered.  In  Ref.  5,  the  differences  between  SM-F  (which  is  close  to  HMS-F) 
and  SL  were  traced  to  the  switching  function  parameters  vi,  y2  and  v3  in  V-.. 

The  values  used  by  SM  (and  by  us)  are  a  factor  of  two  larger  than  those  used 
by  SL.  This  causes  V3  to  be  cut-off  more  rapidly  in  displacements  from 
equilibrium  for  the  MBPT-F  and  HMS-F  potentials  than  SL.  This  explains  why 
the  different  quartic  potentials  used  by  MBPT-F  and  HMS-F  have  so  little 
influence  on  the  shapes  of  tne  curves  in  Fig.  6.  Since  SL 1 s  curves  aareed 
well  with  HML's  ab  initio  results  in  this  region,  it  seems  likely  that  the 
SL  curve  is  more  accurate  than  MBPT-F  and  HMS-F  in  this  region.  This  differ¬ 
ence  does  not  have  a  strong  effect  on  the  thermal  rate  constants,  however.  At 
300  K,  SL  found  that  thermal  rate  constants  for  their  surface  were  about  a 
factor  of  2  lower  than  SM-F.  Since  both  results  were  within  experimental 
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uncertainties,  no  definitive  statement  concerning  the  relative  accuracies  of 
the  two  surfaces  can  be  made  based  on  this  comparison.  The  larger  value  of 
the  SM-F  rate  constant  is  consistent  with  the  longer  range  of  the  attractive 
part  of  that  potential  surface.  This  is  very  likely  to  be  true  with  the 
MBPT-F  surface  as  well. 

It  is  also  of  interest  to  examine  the  behavior  of  the  potentials 
for  linear  OHH  configurations.  Here,  in  agreement  with  Ref.  5,  we  find  that 
the  SM-F,  HMS-F  and  MBPT-F  surfaces  exhibit  wells  while  the  SL  surface  has  a 
small  barrier  followed  by  a  monotonic  decrease  in  energy  along  the  reaction 

0  o 

Dath  going  to  OH  +  H.  The  MBPT-F  well  occurs  at  Ri  =  1.054A,  R3  =  1.005A, 
and  has  an  energy  of  -5.961  eV  relative  to  0(3P)  +  2H ( 2 S )  (1.329  eV  below  the 

O 

0H(- ~  )  +  H  arrangement  channel  energy).  The  HMS-F  well  occurs  at  R,  =  1.C29A, 

R3  =  1.016A,  with  an  energy  of  -6.076  eV  (1.444  eV  below  OH  +  H).  The  SM-F 

e  °  e 

well  occurs  at  Rj  =  1.02A,  R3  =  1.02A,  and  is  1.47  eV  below  OH  +  H' .  The  SL 

o  o 

barrier  is  at  Rj  =  1.18A,  R3  =  0.76A,  and  has  an  energy  1.3  kcal/mole  above 
0(iQ)  +  H2  (2.23  eV  above  OH  +  H).  The  wells  observed  for  the  MBPT-F,  HMS-? 
and  SM-F  are  all  well  above  H20  equilibrium,  but  they  are  Drobably  artifacts, 
as  no  ab  initio  or  semi-emDirical  evidence  supports  their  existance.3  The 
MBPT-F  well  is  only  slightly  smaller  than  HMS-F,  which  indicates  that  neither 
the  additional  quartic  terms  present  in  the  MBPT-0  potential  nor  the  more 
accurate  MBPT-Q  force  constants  has  much  influence  on  this  artifact. 

IV.  CONCLUSION 

In  this  paper,  the  accurate  MBPT  quartic  force  field  of  Bartlett, 
Shavitt  and  Purvis  has  been  fit  to  an  analytical  function  using  a  method 
developed  by  Sorbie  and  Murrell.  Properties  of  this  surface  have  been 
analyzed  and  compared  with  those  of  previously  determined  surfaces  (including 


a  reparametrized  SM  fit  to  the  HMS  force  field).  Our  analysis  indicates  that 
this  fitted  MBPT  surface  describes  the  properties  of  the  H.O  molecule 
ext  emely  well.  In  addition  to  reproducing  the  MBPT  quartic  parameters  at 
the  equilibrium  geometry,  this  surface  reproduces  the  known  energy  associated 
with  straightening  the  H20  molecule  better  than  any  other  global  surface.  Far 
away  from  equilibrium  however,  this  surface  is  much  less  accurate.  For  the 
0(0)  +  H„  arrangement  channel,  the  surface  shows  a  much  longer  range  attractiv 
ootential  than  has  been  observed  in  ab  initio  calculations,  although  the  error 
in  thermal  rate  constants  introduced  by  this  feature  is  apparently  inside 
experimental  uncertainties.  In  addition,  the  0  +  H~  linear  geometry  exhibits 
a  sourious  minimum.  These  errors  all  arise  for  geometries  where  the  3-body 
part  of  the  potential  has  largely  been  damped  out,  so  that  the  more  accurate 
quartic  force  field  used  to  represent  H^O  in  this  surface  has  no  corrective 
influence.  Evidently  to  improve  upon  the  surface  in  these  regions  it  will  be 
necessary  to  use  higher  terms  than  quartic  in  the  H20  force  field,  or  an 
improved  choice  of  the  damping  coefficients  Y:,  y2  and  v3,  or  the  explicit 
incorporation  of  ab  initio  points  far  removed  from  H„0  equilibrium  in  the 
fitting  algorithm. 

Despite  the  apparent  inadequacies  of  this  surface  in  the  O('-D)  + 
regions,  it  does  provide  an  excellent  representation  of  the  H„0  potential 
close  to  equilibrium  (probably  for  energies  as  high  as  3-4  eV  above  the  H.O 

i. 

minimum).  As  such,  this  potential  should  be  useful  for  spectroscopic  studies, 
and  as  the  intramolecular  potential  in  nonreactive  collision  problems  involving 
H20  as  one  partner.  Indeed  for  both  such  studies,  the  MBPT-F  surface  is  much 
to  be  preferred  over  its  quartic  counteroart,  for  important  deviations 
between  the  quartic  and  full  surfaces  can  occur  for  bending  displacements  only 
0.7  eV  above  H20  equilibrium. 
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The  comparisons  of  the  MBPT-F  surface  with  others  indicated  that 
it  is  very  similar  both  qualitatively  and  quantitatively  to  HMS-F  (the 
latter  being  close  to  but  not  quite  the  same  as  SM-F).  Near  H20  equilibrium, 
the  MBPT-F  is  qualitatively  similar  to  but  quantitatively  different  from  the 
SL  surface,  with  SL  the  less  accurate,  while  far  from  equilibrium  (in  the 
0  +  H2  region),  it  appears  that  SL  is  the  more  accurate  surface.  Of  course, 
the  exact  surface  characteristics  in  the  far  from  equilibrium  configurations 
are  still  not  known  very  auantitati vely. 
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TABLE  I.  POLYNOMIAL  COEFFICIENTS  IN  3  BODY  FITTING  TERM 


Terrr 

Pin[S5<Si(cusp)] 

Pout[S,>S,(cusp)] 

5M 

HMS 

MBPT 

SM 

HMS 

MBPT 

1 

1 

1 

1 

1 

1 

1 

S,(S2) 

1.5421 

1.5402 

1.5207 

1.8750 

1.8676 

1 .8559 

S3 

4.6539 

4.6551 

4.6455 

-0.3691 

-0.3582 

-0.35631 

S;2(S22) 

1.5720 

1.5748 

1.2746 

3.1767 

3.1773 

3.0102 

S32 

-4.4355 

-4.4369 

-4.8014 

-2.4168 

-2.4096 

-2 . 6085 

S  ;  S  3  (  S  -  S  3  ) 

17.9914 

17.995 

18.049 

3.1803 

3.1849 

3.2314 

s  :s- 

-1.9761 

-1.9861 

-1.7863 

1 .2577 

1.2402 

1.3496 

S-.?(S23) 

6.5279 

6.0132 

4.3146 

7.1737 

6.8763 

5.9441 

S33 

-7.5781 

-7.5742 

-9.3266 

-5.2200 

-3.2016 

-6.1 567 

S 1 S2  2  ( S  2  S  3  *■ ) 

4.9240 

5.0266 

10.088 

8.1229 

8.1405 

1C. 893 

S  3  S  32  (  S2  s  3-  ) 

11.3559 

11.391 

14.816 

6.2908 

6.2925 

8.1619 

S:2S3(S23S3) 

22.3568 

22.279 

21 .204 

-2.2456 

-2.2399 

-2.7979 

S  3  S2  S  3 

23.1978 

23.089 

12.017 

-1.7899 

-1 .8101 

-7.8075 

S:“(S2-. 

0.6787 

-0.5793 

-3.6377 

0.3677 

-0.3056 

-1.9703 

S  3  “ 

-11.4851 

-11.478 

-14.031 

-6.7255 

-6.7552 

-8.1599 

S:3S;(S;3S:) 

1.5255 

0.5220 

-2.5556 

0.8266 

0.3061 

-1.3837 

S :  - S2  2 

-5.7924 

-5.9434 

4.2511 

-3.1387 

-3.2029 

2.3532 

S’  3  S  3  (  S3  3S  :  ) 

59.6057 

59.011 

63.497 

9.1025 

8.8604 

11.376 

S:2S32  (S22S32) 

-22.2496 

-22.368 

-21.790 

-15.1988 

-15.311 

-15.002 

s  1  s  3 5  ( S2  S  0 3 ) 

20.8998 

20.982 

23.152 

8.0856 

8.1767 

9.3882 

S;S2S3- 

-18.3900 

-18.712 

-20.269 

-13.1075 

-13.316 

-14.186 

S:2S2S3(S:S22S3) 

93.5650 

93.800 

86.304 

15.9060 

16.160 

12.198 

Ain(eV) 

-0.9418 

-0.94159 

-0.94354 

A0Ut(eV) 

-1.7381 

-1.7346 

-1.7344 

S3(cusp)(A) 

0.1568 

0.15677 

0.15638 

Y1  »  Y2  ( ^  ^  ) 

4.5348 

4.5348 

4.5348 

*3  (A-1) 

2.0 

2.0 

2.0 

Vh<!> 

0.9572 

0.9572 

0.95680 

RHp(A) 

1.5139 

1.5139 

1.51429 

■> 

s 

-- 

-- 

25.0 
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TABLE  II.  QUADRATIC,  CUBIC  AND  QUARTIC  FORCE  FIELD 
PARAMETERS  FOR  H20  SURFACES 


HMS-Q 

HMS-F 

SM-F 

MBPT-Q 

MBPT-F 

SL 

wv  1 

3831.5 

3831.4 

3832.0 

3865.0 

3864.9 

3530.2 

'-JO 

1648.8 

1648.8 

1648.6 

1687.4 

1687.4 

1568.2 

‘-'i 

3942.2 

3942.1 

3942.6 

3975.0 

3975.0 

4039.2 

-302.5 

-302.5 

-304.9 

-304.2 

-304.2 

-270.7 

k2o  2 

63.6 

63.6 

63.6 

42.9 

42.9 

10.1 

k:;2 

167.4 

167.4 

168.2 

148.6 

148.6 

235.7 

k::  l 

-53.1 

-53.1 

-53.7 

-61.8 

-61  .8 

-49.0 

k ,  :  3 

-927.8 

-927.8 

-935.5 

-914.1 

-914.1 

-840.9 

k;  3  3 

-138.8 

-138.8 

-139.2 

-111.7 

-111.6 

-35.9 

k  ;  ;  ;  : 

31.9 

31.8 

31 .8 

31.5 

31  5 

34.6 

k2Z22 

2.1 

2.1 

2.1 

-2.6 

-2.6 

25.2 

k  3  3  3  3 

35.4 

35.4 

35.3 

32.0 

32.0 

12.0 

k  •  ?  2 

-85.6 

-85.4 

-86.2 

-75.1 

-74.8 

-65.3 

k  i :  3  3 

201.3 

201.4 

201.3 

190.3 

190.4 

152.5 

k2 2  3  3 

-101.1 

-100.9 

-101.5 

-91.7 

-91.4 

-98.0 

TABLE  III.  EQUILIBRIUM  PROPERTIES  OF  H20  MOLECULE 


Surface 

VH20^ 

b  ° 

R?(A) 

sS<‘ 

SM-F,HMS-F, Experiment 

-10.0705 

0.9572 

104. 

MBPT-F 

-10.0705 

0.9568 

104, 

SL 

-8.95 

0.9867 

103, 

WMF 

-9.16 

0.9778 

92 

Tully 

-9.37 

0.979 

100, 

(a) 

(b) 

(c) 


Potential  at  H,0  equilibrium  geometry  [relative  to  0(;P)  +  H(-S) 

-  +  tj  ( s ) 

OH  equilibrium  position 
H20  angle 


TABLE  IV.  PROPERTIES  OF  LINEAR  -OH  SADDLE  POINT 


Surface 

VeV> 

b 

LV(eV) 

C  o 

R(A) 

d 

ui  (cm-"1 ) 

e 

(cm"1 ) 

f 

w 3 ( cm” 1 ) 

HMS-C 

-8.9G2 

1.168 

0.9506 

3577 

1 407  i 

4064 

MBPT-F 

-8.774 

1.296 

0.9480 

3595 

1  502  i 

4005 

SL 

-7.083 

1.708 

0.9469 

4399 

1 867  i 

4618 

(a)  Potential  at  saddle  point 

relative  to  0( 3P)  +  H(  :S 

)  +  H  ( 1 S ) 

(b)  Difference  between  saddle  point  energy  and  corresponding  H?0  eauilib 
rium  energy 

(c)  OH  distance  at  saddle  point 
(c)  Symmetric  stretch  frequency 
;’e;  Bend  frequency 

(f:  Asymmetric  stretch  frequency 


FIGURE  CAPTIONS 


FIGURE  I.  MBPT-Q,  MBPT-F  and  SL  potentials  (in  eV)  versus  HOH  bend  angle  2 
for  Rls  R2  fixed  at  their  equilibrium  values. 

FIGURE  2.  Equipotential  contours  of  Vh2o  (MBPT-F)  versus  Ri,  R;  for  2 

fixed  at  its  eaui librium  value.  Contours  are  in  1  eV  increments 
starting  with  the  lowest  at  -10  eV  relative  to  0(3P)  +  2H(2S). 

FIGURE  3.  Contours  of  Vh,q  (MBPT-F)  analogous  to  Fig.  2,  but  for  2  =  130° 
(linear  HOH). 

FIGURE  4,  MBPT-Q,  MBPT-F  and  SL  potentials  (in  eV)  versus  Ri  for  R,  =  R2 
and  2'  at  its  equilibrium  value. 

FIGURE  5.  Contours  of  Vh~0  (MBPT-F)  versus  X  and  Rs  for  Derpendi cul ar  0  A  H 
geometries.  The  distance  X  is  the  0  to  center  os  mass  of  H2  dis¬ 
tance.  Contours  chosen  are  the  same  as  in  Fig.  2. 

FIGURE  6.  MBPT-F,  HMS-F  and  SL  DOtentials  versus  X  for  perpendicular  0  +  H; 
and  R3  equal  to  the  H2  equilibrium  distance. 

FIGURE  7.  MBPT-F,  HMS-F  and  SL  potentials  versus  X  as  in  Fig.  6,  but  for  R-. 
equal  to  the  H20  equilibrium  H2  distance. 
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Abstract 

/V\.'WW\A> 


We  calculate  the  adiabatic  potential  energy  curves  and  nonadiabatic 
first-derivative  couplings  for  the  X,  A,  and  C  states  of  KH  by  an  Ab  initio 
one-electron  pseudopotential  formalism.  The  splitting  of  the  X  and  A  curves 
at  the  avoided  crossing  is  in  good  agreement  with  experiment.  The  ab  initio 
results  are  used  to  calculate  the  electronically  inelastic  transition  proba¬ 
bilities  and  cross  sections  for  K  +  H  collisions  at  low  energies  by  R  matrix 
propagation  in  the  adiabatic  representation  with  exponential  sector  trans¬ 
formations.  Oince  this  method  has  never  been  applied  before,  we  made  an 
extensive  study  of  its  convergence  properties  and  efficiency.  We  found 
it  to  be  a  convenient,  accurate,  and  efficient  method.  The  cross  sections 
are  changed  by  about  a  factor  of  two  when  the  potential  curves  are  changed 
by  a  different  treatment  of  the  KH+  core,  but  only  by  about  1%  when  the 

assumptions  about  the  nonadiabatic  second-derivative  coupling  terms  are 

2  2 

altered.  Our  estimate  of  the  4  P  -+■  4  S  quenching  cross  section  at  0.022  eV 

— 4  2 

relative  translational  energy  is  2-4  x  10  a q.  This  increases  to 

-4  2 

8-10  x  10  Sq  by  1.1  eV.  The  emphasis  in  this  article  is  on  testing  and 
evaluating  the  new  method  for  solving  the  scattering  problem  rather  than 
on  the  cross  sections  themselves. 


I.  INTRODUCTION 

The  standard  quantum  mechanical  treatment  for  low-energy  atomic  and 

molecular  collisions  is  the  close  coupling  method.''"  When  more  than  one 

electronic  state  must  be  considered,  one  can  use  an  adiabatic  or  a  diabatic 

representation  for  the  electronic  wavefunction.  The  diabatic  representation 

has  the  mathematical  convenience  of  no  derivative  coupling  operators,  but 

2  3 

it  is  not  unique.  ’  One  way  to  specify  it  completely  is  to  define  it  by 
a  transformation  from  a  finite  number  of  adiabatic  states,  where  the  trans¬ 
formation  is  defined  by  requiring  the  first-derivative  coupling  to  vanish 
in  the  finite  manifold.  It  is  also  possible  to  solve  the  coupling  equations 
directly  in  the  electronically  adiabatic  representation,  including  the 

derivative  coupling.  A  new  method  for  doing  this  has  been  proposed  by 

3 

two  of  the  authors,  and  it  is  applied  here  for  the  first  time.  The  method 
involves  R  matrix  propagation  and  requires  as  input  only  the  adiabatic 

4 

potential  curves  and  first-derivative  coupling  matrix  elements  obtainable 
from  standard  electronic  structure  calculations.  For  the  present  appli¬ 
cation  we  consider  collisions  of  K  with  H,"*  we  consider  only  radial 
coupling  between  '"E+  states,  and  we  obtain  the  adiabatic  potential  curves 
and  first-derivative  coupling  matrices  by  ab  initio  methods. 

Section  II  presents  the  coupled-channels  scattering  equations  and  the 

details  of  how  we  solve  them.  This  section  also  compares  the  new  method  to 

12 

the  method  Johnson  and  Levine  “  proposed  for  this  problem  and  to  a  method 
13 

one  of  us  and  Wyatt  have  applied  to  solve  reactive  scattering  problems 
in  vibrat iona Lly  and  rotationally  adiabatic  representations.  Section  III 
gives  details  and  results  of  the  electronic  structure  calculations  performed 
to  -enerate  the  input  to  the  scattering  equations.  Section  IV  presents  the 
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details  of  the  scattering  calculations,  section  V  presents  results,  and 
section  VI  is  discussion.  The  emphasis  in  the  present  paper  is  on  evaluating 
the  new  method  for  solving  the  scattering  equations  rather  than  on  the 
cross  sections  themselves. 


-3- 


II.  THEORY 


A.  New  method 

Consider  the  coupled  I  states  of  a  diatomic  system  with  nuclear  masses 
m^  and  nig  and  N^  electrons  each  of  mass  m^.  The  total  Hamiltonian  is  given 
by 


H 


+  H  (R) 
e 


2(mA+V 


(1) 


(2) 


where  R  is  the  internuclear  distance,  r.  is  the  vector  from  the  center 

of  mass  of  the  nuclei  to  electron  i  in  the  bodv-fixed  frame,  L  is  the 

angular  momentum  operator  of  the  relative  motion  of  the  nuclei,  y  is 

AB 

the  reduced  mass  for  this  motion,  and  H  is  the  "electronic  Hamiltonian" 

e 

+  2  2 

H  (R)  =  -  ~z  E  v  +  V (x,R) 
e  zm  .  r . 

e  i  i 

Ne 

where  we  have  denoted  x  as  the  collection  of  electronic  coordinates  <r.;.  ,, 

’  l  1=1’ 

and  V(x,R)  includes  all  the  pairwise  coulomb  interactions  between  nuclei 
and  elec-trons. 

The  radial  coordinate  is  subdivided  into  sectors  numbered  (i),  and 
within  each  sector  the  total  wavefunction  Y^(x,R),  where  q„  denotes’  the 
initial  conditions,  is  expanded  in  an  orthogonal,  approximately  adiabatic 
basis  ti  (jc;i)  which  is  independent  of  internuclear  coordinate: 


=  [pa(^c,i)  ]T  xa(R;i) 


(3) 


where  ?  is  a  row  vector  of  elements  ?  ,5  is  a  column  vector,  each  column 

a  '  qo  * 

of  <  is  a  different  linearly  independent  solution,  and  each  row  of  \ 


corresponds  to  one  of  the  channels. 
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The  approximately  adiabatic  basis  states  in  sector  (i)  are  chosen  to 
diagonalize  the  molecular  electronic  Hamiltonian  in  the  body-fixed  frame 
at  the  center  of  the  sector,  i.e.. 


/d?»‘(?il>“X>  <fS5=1>  •  Vaeo(1> 


q1q  q 


(4) 


Coupled  radial  equations  for  each  sector  (i)  are  obtained  by  substituting 

eqs.  (l)-(3)  into  the  time-independent  Schrodinger  equation  for  total  energy 

E,  and  closing  on  the  left  with  an  adiabatic  basis  function  4>a  (x,i).  Neg- 

qi  % 

lecting  the  mass  polarization  term  [last  term  of  eq.  (1)]  and  Z-H  coupling, 
the  coupled  equations  in  sector  (i)  are 
>2  2 

{-  —  [ — y - y-*- 1  +  t  (x)  -  E }  X-  _  CR;i) 

2uAB  dR2  R2  qi  Vo 


-vv(R;i)H(R;i)  =  0 


(5) 


where 


V”  (R;  i)  =  /dx  ia  (x;  i)  H  (R)  i>a(x;i) 

'  '  n  rV,  q  q 


qLq 


ql  % 


(i)  3 


qiq 


(6) 


The  matrix  elements  V  (R,'i)  are  defined  for  all  values  of  R  in  sector  (i) 

qiq 


using  basis  functions  that  would  usually  be  used  only  for  calculations  at 
the  center  of  the  sector;  these  matrix  elements  are  zero  at  the  center  R* 
of  each  sector.  Also,  because  the  basis  functions  fa(;<;i)  are  independent 
of  R  within  a  sector,  the  coupled  radial  equations  in  each  sector  contain 
no  derivative  coupling  terms.  A  more  convenient  form  of  eq.  (5)  is 


:<a  ( R ;  i)  =  ' 2  (R ;  i)  ,a(R;i)  (7) 

dR“ 


where 
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A“(R.  i)  =  iVd(R;i) 

tr  J 


EaCi)  + 

2'^abr' 


-  E  ]  I  } 


and  Che  diagonal  matrix  E  (i)  is  given  by 


E3  (i)  =  ca(i)  6 

qq,  q  qq. 


(8) 


(9) 


Using  standard  numerical  techniques,  the  radial  wavefunctions  x3(R;i)  can 

a 

be  propagated  from  the  left  side  to  the  right  side  of  sector  (i).  Continuous 
solutions  with  continuous  first  derivatives  are  obtained  for  multi-sector 
regions  by  imposing  the  following  matching  conditions  at  each  sector  boundary 


*(1-1) 

’q 


,  „i-l, 

(^’RR  } 


(10) 


dR  q  C$,RR  ]  dR  ‘q  AV 


(11) 


,i,A 


where  R^(R^)  is  the  value  of  R  at  the  right  (left)  boundary  of  sector  (i), 
Substituting  eq.  (3)  into  eqs.  (10)  and  (11)  and  rearranging  yields 


d  a,  i-l  .  T  d  a,  i  .. 

dR  i  (RR  ’1_1)  -v(l  1;l)  dR  'X  (RL;l) 


(12) 


xa(RR_1;i-D  =  T( i-i;i)  xa(R7Si) 

A  K  -  /v  L 


(13) 


where 


T 

qq 


l 


(i-l; i) 


/dx  *a"(x;i-D 


ia  (x;  i) 
ql  J 


(14) 


Equations  (12)  and  (13)  provide  the  initial  conditions  for  propagation 

across  sector  (i)  given  the  solution  in  sector  (i-l).  Thus  a  continuous 

solution  can  be  constructed  over  the  entire  scattering  region.  The  equations 

presented  so  far  are  essentially  t he  same  as  those  in  the  method  proposed 

12 

by  Johnson  and  Levine.  “  However,  instead  of  requiring  an  explicit  evaluation 


of  the  inconvenient  overlap-type  integration  in  (14),  we  relate  the  trans¬ 


formation  matrices  T(i-l,i)  to  the  standard  nonadiabatic  derivative  coupling 
matrices. 

To  accomplish  this,  we  define 


M  (R„,y)  = 
sq  C  3 


«t> 


{i'Kc + 


S  "V. 


(15) 


where  the  matrix  element  indicates  integration  over  the  electronic  coordi¬ 
nates.  Then 


•  h-  +h-+i 

T  ( i ,  i+1 )  =  M(Rg,  1'2-  ) 


(16) 


where  tu  is  the  length  of  sector  (i).  Differentiating  (15),  we  obtain 


5y 


M 

sq 


,y) 


<Oa(x,Rb!-r^| 

S  -j  C  ay 


a ,  i 

(x»Rr 

q  v  C 


y)> 


,y) 


+  y) 


where  F  (r)  is  the  nonadiabatic  derivative  coupling  matrix  defined  by 


(17) 

(13) 


(19) 


In  obtaining  eqs.  (18)  and  (17)  we  assumed  that  4-a(x,i)  form  a  complete  set 
of  states.  When  this  approximation  is  used  in  coupled-channel  calculations 
employing  a  truncated  set  of  states,  its  validity  can  be  '  .sured  by  obtaining 
converged  results  with  respect  to  increasing  the  basis  size.  The  effects 
of  this  approximation  for  finite  basis  sets  are  discussed  further  in  the 
Appendix.  Equation  (18)  can  be  solved  numerically  for  M(R*,y)  using  the 
Magnus  method, 


i 
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M(R^,v)  =  M(Rj,0)exp  [N(R^,y)] 

where 

.  2u  y 

N(Rp,y)=  -  -p  /  dy’  Fa(R^  +  y') 
^  C  ti2  0  J  C 


(20) 


and 


+  2(-4fr  /'dy'  /  dy" [Fa (R^  +  y’),Fa(Ra  +  y")  ]  + 
Vo  o  J  J  L 


(21) 


«  (RC'0)  =  l 


(22) 


*M  "  «RC 


Ri+l> 


(23) 


and  expand  Che  coupling  matrix  in  a  Taylor  series  around  this  point 

.  dFa 

f (Rj  +  y)  -  f(*i>  +  (y  +  Rj  -  *J)  ir(i* 


(24) 


hi  +  hi+l 

Substituting  (24)  into  (21)  and  retaining  terms  through  order  - - -  yields 


•  h .  +  h 

i  i  l+l 


2(RC’  2 


)  = 


2uD  h .  +  h 

AB  ,  l  i+l  .  _a l 


h.  +  h  . 

l  i+l s  3. 


(  2  >  r(v  +  °u-— — — )  i 


(25) 


and,  by  (16)  and  (20), 


2u  h .  +  h  .  ^ .  +  h  .  - 

t  , .  ,  AB  i  i+l  „a,„iN,  .  l  i+1n3, 

|(i,  i+l)  =  exp  [ - 2 - 2 -  l  °K - 3 - )  ] 


(26) 


This  completes  the  derivation  of  a  way  to  perform  the  sector-boundary  matching 
of  (12)  and  (13)  by  using  the  standard  functions  jfa(R)  rather  than  the  non¬ 
standard  overlaps  of  (14). 

The  essential  step  in  tiie  new  method  is  eq.  (26)  for  the  sector  trans¬ 
formation  matrix.  We  note  that  this  expression  is  identical  to  the  sector 


w 
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transformation  matrix  given  in  eq.  (12)  of  reference  13.  The  derivation 
and  the  context  are  different,  however.  In  reference  13  the  problem  con¬ 
sidered  was  reactive  scattering  in  a  vibrationally-rotationally  adiabatic 
basis,  and  the  propagation  of  the  radial  wave  functions  of  (7)  from  to 
Rd  was  accomplished  by  the  first  order  Magnus  method . ^ ^  In  the  present 
case  the  internal  degrees  of  freedom  being  considered  are  the  electronic 

ones,  and  we  will  accomplish  the  propagation  across  a  sector  by  the  R  matrix 
3  16-18 

propagation  method.  ’  Comparison  of  the  approaches  illustrates  that 

eq.  (26)  provides  a  general  sector  transformation  matrix  for  solving  scat¬ 
tering  problems  in  adiabatic  bases.  For  problems  involving  both  electronic 
and  vibrational-rotational  degrees  of  freedom,  such  as  electronically  inelastic 
atom-molecule  or  molecule-molecule  scattering,  it  is  possible  either  (i)  to 
use  (26)  or  eq.  (12)  of  reference  13  to  treat  all  degrees  of  freedom  in 
an  adiabatic  representation  or  (ii)  to  treat  electronic  degrees  of  freedom 

in  an  adiabatic  representation  and  to  treat  other  degrees  of  freedom  by 

19 

the  diabatic  representation  that  is  more  standard  for  vibrational-rotational 
degrees  of  freedom. 


B.  Standard  methods 

To  compare  with  the  adiabatic-at-the-center-of-a-sector  method  given 

above,  we  present  here  the  equations  for  the  standard  adiabatic  and  diabatic 
.  ,  20 

propagation  methods. 


Using  a  continuous  adiabatic  basis  J>  (x,R)  the  coupled  equations  analogous 

n,  % 

to  eq.  (5)  are 


*2  H2 

n  f  d _ 

“  “AB  dR“ 


-j1-  +  — ^r~]  +  fa(R)  +  2Fa (R)  ~  +  G3(R)  >x*(R)  =  0 
r2  ti  dR  " 
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where 


E'l  (R)  =  ef(R)  5 


2 1  a 
3R  ql 


and  F  (R)  is  defined  in  (19).  Equation  (27)  contains  the  radial-derivative 
couplings  Fa(R)  and  Ga(R),  but  the  angular-derivative  couplings  do  not  appear 
because  we  restricted  ourselves  to  Z  states  at  the  beginning  of  section  II. A. 
Note  that  x  (R;i)  and  X&(R)  become  identical  in  the  limit  of  small  sector 
sizes.  The  first-derivative  coupling  term  can  be  eliminated  from  (27)  by 
transforming  to  an  orthogonal  diabatic  basis  such  that  the  transformation 
matrix  ( R)  obeys  the  following  equation 

Ad 

The  transformed  coupled  equations  are  then  given  by 

t-  ~  •t(*'>21)  +  — r-1  +  %d(R)  +  £d(R>}*dw  = 0  (3d 

^AB  dR  R  h  ^ 


where 


(R)  =  U(R)  Ea(R)  UT(R) 


Xd(R)  =  U(R)  Xa(R) 


and  the  second-derivative  coupling  in  the  new  basis  is 

Gd(R)  =  U(R)  (Ga(R)  -  [Fa(R)]2}UT(R) 

n 

It  can  be  shown  that  in  the  limit  of  a  complete  basis 


S'1(R>  =  -IT 


dFa  2u 


f^[Fa(R)]2 


and  thus  G  (R)  vanishes  in  that  limit.  The  diabatic  basis  defined  by  equations 
(30)  and  (32)  is  a  P-diabatic  basis  in  the  terminology  used  in  references  2 


and  3 . 
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III.  SYSTEM  AND  METHODS 

The  system  studied  for  a  test  case  is  K  +  H.  The  adiabatic  potential 

11  21 

curves  were  calculated  by  a  one-electron  model  ’  for  alkali  hydrides 
involving  effective  core  potentials  to  represent  K+  and  H.  This  method 
is  described  in  reference  11,  and  the  effective  core  potentials  and  orbital 
basis  set  used  are  also  given  there.  The  three  lowest-energy  adiabatic 
states  4>j(x;R)  of  E+  symmetry  were  calculated  by  diagonalizing  the  one- 
electron  Hamiltonian  in  the  orbital  basis.  The  potential  curves  for  these 
states  were  calculated  by  adding  the  energy  of  the  KH+  core  to  the  one- 
electron  eigenvalues.  The  KH+  core  energy  is  approximated  by  a  full  KH+ 
calculation  or  by  a  calculation  on  KH+  employing  only  a  single  H  Is  basis 
function.11  Following  reference  11,  we  call  these  two  choices  methods 
21  and  2H,  respectively.  As  discussed  there,  method  2H  is  expected  to 
be  the  more  accurate  one.  The  derivative  coupling  matrix  Fa(R)  was  then 
calculated  from  the  wave  functions  by  using 

<<?>j(£;R)  =  ?  Mjk(R,8)  (36) 

since  <l!>jl*jc>  "0.  A  value  of  0.001  a^  was  used  for  6.  For  calculating 

the  overlap  integral  M.,  (R,  6)  ,  defined  by  (15),  the  K  nucleus  was  fixed 

J  K 

and  the  H  was  moved  by  the  amount  <5.  This  corresponds  to  placing  the  origin 

of  the  electronic  coordinate  system  at  the  K  nucleus,  and  it  makes  the 

3*3  submatrix  of  £a(R)  for  the  states  considered  here  tend  to  a  null  matrix 

at  R  =  »,  which  makes  it  straightforward  to  impose  scattering  boundary  conditions. 

3  22 

Other  possible  choices  of  electronic  origin  are  discussed  elsewhere.  ’ 

The  calculated  adiabatic  potential  curves  ea(R)  are  shown  in  Figure  1, 
and  the  calculated  first-derivative  coupling  matrices  are  shown  in  Figure  2. 

In  Figure  2  we  use  the  notation 
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Fa(R)  =  -(^i2/2u  )  fa(R)  (37) 

x  AB 

Transforming  away  the  first-derivative  coupling  by  the  3  *  3  matrix  U(R) 
yields  the  diabatic  Hamiltonian  matrix  H^(R)  whose  elements  are  illustrated 
in  Figures  3  and  4. 

The  results  are  given  in  hartree  atomic  units:  1  a.u.  energy  =  1  hartree  = 

Q 

1  E,  =  27.212  eV,  and  1  a.u.  length  =  1  bohr  =  1  a^  =  0.52918  A. 
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IV.  SCATTERING  CALCULATIONS 

We  will  compare  several  different  methods  for  solving  the  scattering 
problem. 

A.  First  order  Magnus  approximation  in  the  adiabatic  representation. 
In  this  method  we  rewrote  (27)  as 


dY 

■v 


dF  "  £<R>  l  <R> 


(38) 


where 


f(R)  = 


I 

'V/ 


2fa(R) 


(39) 


(40) 


and 


J8(R)  =  ^[Ea(R)  +  <f(R)  +  (^iii±ll  -  E)  %] 


2*abr' 


(41) 


Equation  (38)  was  integrated  by  the  first  order  Magnus  approximation 


Ya(R  +  h)  =  exp[hA(R  +  i5h)]^a(R)  (42) 

The  exponential  was  evaluated  by  a  power  series,  retaining  terms  through  h  . 
Checks  showed  that  the  same  results  were  obtained  by  retaining  terms  through  h\ 

The  scattering  matrix  was  evaluated  by  applying  boundary  conditions  to 
a  23 

X  (R)  in  the  usual  way.  We  employed  a  fixed  stepsize  h  and  decreased  it 
% 

till  convergence  was  obtained  for  the  absolute  squares  of  scattering  matrix 


elements. 


with  power  series  evaluation  of  the  exponentials.  In  this  method  we  rewrote 


(31)  as 


“d F  =  S(R)  l  (R) 


where 


Yd(R)  = 


xd(R) 

% 


Dd(R) 


AB  r„d, 


2 

fW  =  -f-  [H“ (R)  +  jg“(R)  +  (  —  t-p-  -  E)  I  ] 


2\br 


S  S' 

Starting  with  the  same  input,  E  (R)  and  F  (R) ,  as  for  method  A,  we  integrated 
(30)  simultaneously  with  (43)  so  that  we  could  calculate  B(R)  from  (32),  (45), 
and  (46)  as  we  needed  it.  We  assumed  Gd  =  0  as  discussed  after  (34).  The 
first  order  Magnus  approximations  are 

UT(R  +  4h)  =  exp[-hfa(R)]  UT(R- ^jh)  (47) 


Yd(R  +  h)  =  exp[hB(R  + 'rh)  ]Yd(R)  (43; 

In  (47)  and  (43)  the  exponentials  were  evaluated  by  power  series  through  h2 . 


The  boundary  conditions  and  stepsize  were  handled  the  same  as  in  method  A. 
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C.  Firsc  order  Magnus  approximation  in  the  diabatic  representation 
with  analytic  evaluation  of  the  exponentials.  This  method  is  the  same  as 
method  B  except  for  the  evaluation  of  the  exponentials.  The  exponential 
in  (48)  was  evaluated  analytically  in  terms  of  the  eigenvalues  and  eigen¬ 
vectors  of  D^(R)  as  explained  by  Light. ^  Since  we  still  assume  G^(R)  =  0, 
the  eigenvalues  are  the  already-available  ea(i)  and  Che  eigenvectors  are 
the  columns  of  U(R  ).  The  latter  are  obtained  from  (47),  but  in  this  method 

*  c 

the  exponential  in  (47)  was  evaluated  analytically.  Since  f  is  skew  sym- 

SL  cL 

metric,  the  exponential  in  (47)  can  be  evaluated  in  terms  of  f^,  f^>  and 

f a  •  This  kind  of  procedure  is  very  efficient  for  2><2  and  3x3  cases,  and 

we  used  it  for  the  calculations  reported  here.  For  matrices  of  order  greater 

than  3,  the  exponential  of  a  skew  symmeteric  matrix  can  be  evaluated  efficiently 

24  25 

by  diagonalizing  the  square  of  the  skew  symmetric  matrix.  ’  The  boundary 
conditions  and  stepsize  were  handled  the  same  as  in  method  A. 

D.  R  matrix  propagation  with  adiabatic  basis  functions  at  sector  centers 

and  with  the  exponential  sector  transformation  matrix.  This  is  the  new 
method  of  reference  3  and  the  present  paper.  The  propagation  across  a 
sector  was  accomplished  by  the  R  matrix  propagation  method, ^  using 
a  modified  version  of  an  R  matrix  propagation  code  that  has  been  discussed 
elsewhere. ^>27  5£nce  Fa(R)  is  skew  symmetric,  the  exponential  sector  trans¬ 
formation  matrix  (26)  was  evaluated  analytically  in  terms  of  Fa?,  and 

cl 

^22,  as  discussed  in  subsection  C  above.  The  method  for  extracting  a  scat- 

27 

tering  matrix  from  the  global  R  matrix  is  explained  elsewhere.  We  used 

26  27  ( ]_ ) 

a  variable  stepsize  algorithm  ’  with  one  stepsize  parameter  s  for 

R  <  25  aQ  and  another  e ^  for  R  >  25  aQ.  was  set  at  a  value  that 

yields  high  accuracy,  and  was  decreased  till  convergence  was  obtained 

for  absolute  squares  of  scattering  matrix  elements. 
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E,  F.  Numerical  integration  in  the  adiabatic  representation.  For  these 

28 

methods  we  applied  a  f ixed-stepsize  Runge-Kutta-Gill  integration  (method  E) 

9  Q 

or  a  5th  order  variable-stepsize  predictor-corrector  algorithm"  (method  F) 
to  integrate  eq.  (38).  The  boundary  conditions  were  handled  the  same  as  in 
method  A. 

G,  H.  Numerical  integration  in  the  diabatic  representation.  Finally 
we  applied  the  Runge-Kutta-Gill  (method  G)  and  5th  order  variable-stepsize 
predictor-corrector  (method  H)  methods  to  simultaneously  integrate  eqs. 

(30)  and  (43) .  The  boundary  conditions  were  handled  the  same  as  in  method  A. 

I.  Initial  values  and  boundary  conditions.  For  all  the  methods  we 
started  the  integration  of  the  radial  wavefunctions  at  small  enough  R  that 
the  results  are  invariant  to  further  decreasing  the  starting  value.  In  par¬ 
ticular  we  started  the  s-wave  solutions  at  1.1  aQ  and  2.1  aQ  for  the  21 
and  2H  potential  curves,  respectively,  for  E  =  0.06  E^.  For  higher  E  we 
started  at  smaller  R,  and  for  higher  l  we  started  at  larger  R,  e.g.,  2.2  aQ 
for  l  =  15  for  the  2H  potential  curves  for  E  =  0.06  E,  .  At  the  starting 
point  the  radial  wavefunctions  were  taken  as  zero  and  the  matrix  of  radial- 
wavefunction  derivatives  was  the  unit  matrix.  This  generates  N  linearly 
independent  solutions  where  N  is  the  number  of  states  retained  in  the  wave- 
function  expansion.  At  R  =  70  a^,  we  transformed  to  the  adiabatic  repre¬ 
sentation  (only  necessary  in  methods  B,  C,  G,  and  H)  and  took  linear  com¬ 
binations  of  the  linearly  independent  solutions  to  obtain  the  correct 
scattering  solutions  satisfying  the  Ricatti-Bessel  boundary  conditions 
for  the  radial  wavefunctions  and  their  derivatives. 

For  methods  B,  C,  G,  and  H  we  need  also  specify  the  initial  values  for 
T 

U  at  the  center  of  the  first  sector.  It  can  be  shown  that  the  scattering 


matrix  obtained  by  the  above  procedures  is  invariant  to  the  starting  value 
T 

of  y  so  we  started  it  as  the  unit  matrix.  This  transformation  generates 

an  arbitrary  linear  combination  (with  R-independent  coefficients)  of  the 

physical  diabatic  states,  where  the  physical  ones  are  the  ones  that  make 
^  T 

P  (R)  diagonal  at  large  R.  Let  C  (R)  denote  the  transformation  matrix 

T 

for  the  physical  diabatic  states .  C  (R)  is  not  needed  for  scattering  cal- 
culations  but  is  required  to  make  Figures  3  and  4.  It  is  generated  by 

CT(R)  =  U(large  R)  L’T(R)  (49) 

r\j  'Xj  'Xi 

where  large  R  =  70  a^  in  practice  and  both  matrices  on  the  right  side  of 

T 

(49)  are  generated  by  solving  (47)  with  U  (R  +  ^h)  =  I. 

a.  start  a- 

In  all  cases  we  checked  that  ending  the  integrations  at  R  =  69  aQ 
would  have  given  the  same  scattering  matrices  (within  0.1%)  as  ending  at 
R  =  70  . 

J.  Stabilizing  transformations.  For  all  the  methods  except  R  matrix 

propagation,  method  D,  it  is  necessary  to  integrate  the  radial  wavefunction 

through  classically  forbidden  regions.  In  such  regions  one  faces  the  well 

known  problem  that  components  of  the  solution  vectors  grow  exponentially 

30-32 

and  can  cause  the  solution  vectors  to  become  linearly  dependent.  This 

problem  is  handled  by  performing  a  Schmidt  orthogonalization  of  the  solution 
matrix,  X3(P0  or  y^(R),  after  a  specified  number  of  integration  steps  have 
been  taken  in  regions  in  which  at  least  one  channel  is  closed.  For  step- 
sizes  small  enough  to  insure  0.01%  accuracy  it  was  found  that  stabilizing 
every  20  steps  was  sufficient  so  that  the  transition  probabilities  were 
invariant  to  6  significant  figures  to  or thogonalizing  even  more  often.  Doublinc 
the  number  of  orthogonalizations  increased  the  computer  time  bv  onlv  10%. 

For  the  larger  stepsizes  which  give  approximately  0.5%  accuracy,  stabilizations 


were  performed  every  step  in  regions  in  which  at  least  one  channel  is  closed. 
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V.  RESULTS  AND  DISCUSSION 

In  all  cases  we  included  three  states  in  the  coupled-channel  calcu¬ 
lations  . 

A.  Comparison  of  methods.  We  made  a  detailed  study  of  computational 

efficiency  for  the  case  of  the  s-wave  probabilities  at  E  =  0.06  E^.  First 

we  performed  some  calculations  with  very  small  stepsizes  (i.e.,  small  sector 

widths)  to  get  the  accurate  transition  probability  P  ^  connecting  the  first 
2  2 

(4  S)  and  second  (4“P)  atomic  states.  Then,  for  the  six  most  efficient 
methods  we  performed  calculations  for  a  fine  grid  of  fixed  stepsizes  (methods 
A,  B,  C,  E,  and  G)  or  fixed  stepsize  parameters  (method  D)  to  find  the 

minimum  computing  time  required  to  achieve  0.5%  accuracy  for  this  probability. 

These  computing  times  are  given  in  Table  I,  where  they  are  expressed  as 
ratios  to  the  computing  time  required  by  method  D.  We  see  that  the  new 
method,  i.e.,  the  variable-stepsize  R  matrix  propagation  method  in  the 
adiabatic  representation  with  the  exponential  sector  transformation  matrix, 
is  the  most  efficient  of  the  seven  methods  tested  in  this  work.  Second 
most  efficient  is  f ixed-stepsize  Magnus  integration  in  the  diabatic  repre¬ 
sentation  with  analytic  exponentiation.  Analytic  exponentiation  was  about 
twice  as  fast  as  using  the  power  series,  even  though  the  power  series  was 
coded  very  efficiently  to  take  advantage  of  the  structure  of  zeroes  in 
the  matrices.  The  Magnus  methods  could  probably  be  made  more  efficient 
by  using  a  suitable  variable-stepsize  algorithm,  but  this  was  not  attempted. 

It  is  well  known  that  the  most  efficient  method  for  one  level  of  accuracy 
is  not  necessarily  the  most  efficient  method  for  other  levels  of  accuracy. 

A  more  detailed  comparison  of  the  two  most  efficient  methods,  on  the  same 
scale  as  used  for  Table  I,  is  given  in  Table  II.  Table  II  shows  that  the 

I 

} 


convergence  of  both  methods  is  smooth  and  that  method  D  is  also  efficient 
for  higher  accuracy,  e.g.,  it  is  3.7  times  more  efficient  than  method  C  for 
0.1%  accuracy  in  • 

It  should  be  clear  that,  although  eight  methods  have  been  compared 
for  the  same  problem,  there  has  been  no  attempt  to  determine  the  absolutely 
fastest  possible  way  to  solve  the  coupled-channels  problem.  The  main  con¬ 
clusion  of  this  section  is  that  the  R  matrix  propagation  method  in  the 
adiabatic  representation  with  the  exponential  sector  transformation  matrix, 
which  is  a  very  convenient  and  stable  method,  is  also  very  efficient. 

B.  Probabilities.  The  s-wave  transition  probabilities  tor  the  two 
sets  of  potential  curves  are  shown  in  Figure  5.  These  probabilities  show 
regular  oscillations  as  functions  of  1/E,  and  the  envelope  of  the  oscil¬ 
lating  probabilities  increases  with  increasing  energy  above  thresw,old. 

The  small  differences  between  the  two  sets  of  potential  curves  change  the 
phase  of  the  oscillations  in  the  transition  probabilities,  but  they  do  not 
make  large  changes  in  the  magnitudes  of  the  envelopes.  For  both  sets  of 

potential  curves  ? ^  is  in  the  range  10-5-10-Zl,  P  is  10~7-10-6,  and 
_  2  _2 

is  10  -  5  *  10  for  most  energies  considered. 

We  performed  some  extra  calculations  to  check  the  assumption  of  a 
complete  set  of  states  for  the  second-derivative  coupling  term.  In  the  standard 
versions  of  methods  A,  E,  and  F,  we  use  equation  (35)  for  Ga .  This  equation 
is  also  used  to  make  0^  vanish  in  methods  B,  C,  G,  and  H,  and,  as  discussed 
in  the  appendix,  it  is  required  to  hold  for  the  equivalence  of  method  D 
•.  >  the  other  methods.  Thus  it  is  interesting  to  test  the  sensitivity  of 
•  r.  ■  ;lt  to  the  treatment  of  the  second-derivative  coupling. 


The  correct  expression  for  a  second-derivative  coupling  matrix  element 


GVR)  =ikFij(a) 


2y 


+ 


AB 


F^CR) 


Fkj(R) 


(50) 


where  the  sum  over  k  should  include  a  complete  set  of  states.  In  (35), 
however,  the  sum  includes  only  the  N  states  retained  in  the  expansion  (3) 
of  the  wavefunction.  To  test  the  importance  of  this  truncation  in  the 
second  term,  we  repeated  the  method-A  calculations  entirely  neglecting 
the  second  term.  The  results  are  given  in  Tables  III  and  IV.  In  27  out 
of  28  cases,  the  difference  of  the  results  is  1%  or  less.  In  the  remaining 
case  the  difference  is  5%.  The  present  test  is  a  very  stringent  one  for 
eq.  (35).  First  of  all,  the  second  derivative  coupling  is  known  to  be 

20  33 

less  important  in  high-energy  cases  where  semiclassical  methods  are  valid,  ~ ’  ’ 

but  the  present  tests  are  low-energy,  highly  quantal  cases.  Second,  the 
truncation  of  the  sum  in  the  second  term  of  (50)  would  be  expected  to  be 
most  valid  for  large  inelastic  probabilities  dominated  by  two  strongly 
coupled  states.  But  the  present  inelastic  probabilities  are  very  small 
and  do  not  correspond  quantitatively  to  a  simple  two-state  avoided  crossing. 

Thus  the  fact  that  the  second  term  of  equation  (35)  has  only  a  small  effect 
in  the  present  cases  is  very  encouraging. 

It  should  be  noted  that  approximations  to  the  skew-symmetric  first 
term  of  equation  (35)  or  (50)  are  dangerous.  If  this  term  is  neglected, 
the  calculated  probabilities  no  lonser  sum  to  unity  or  satisfy  microscopic 
reversibility.  The  second  term  of  (35)  is  symmetric;  thus  its  neglect  does 
not  affect  these  properties. 


C.  Cross  sections.  Cross  sections  were  calculated  at  two  energies, 


E  =  0.06  E,  and  0.10  E.  . 

h  h 


The  first  of  these  is  only  0.022  eV  above  the 
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2 

4“P  threshold  at  E  =  0.059192  E.  ;  thus  it  is  a  total  energy  typical  of 

h 

2 

those  contributing  to  quenching  of  the  4  P  state  under  thermal  conditions. 

2 

The  second  energy  is  1.11  eV  above  the  4  P  threshold.  The  cross  sections 

are  given  in  Table  IV.  At  0.06  E^>  the  potential  curves  obtained  by  method 
2  2 

21  lead  to  4  P 4  S  quenching  cross  section  1.9  times  smaller  than  is 
obtained  by  method  2H.  However  at  0.10  E^,  the  method  21  cross  sections 
are  1.3  times  larger.  Considering  the  small  size  of  the  cross  sections 
and  the  out-of-phase  oscillations  in  the  fixed-X.  inelastic  transition  prob¬ 
abilities,  the  different  values  obtained  for  the  cross  sections  are  not 

2  2 

too  surprising.  The  much  larger  4  P+5  S  cross  section  is  less  sensitive 
to  the  difference  in  the  potential  curves. 

D.  Potential  curves  and  coupling  terms.  Figure  1  compares  the  new 

potential  curves  to  the  experimental  ones.  Method  21  leads  to  more  accurate 

dissociation  energies  but  method  2H  leads  to  more  accurate  repulsive  walls. 

35 

The  C-state  curve  shows  a  shoulder  at  R  =  5-6  a^.  A  similar  feature  was 
predicted  in  reference  7.  Reference  7  shows  that  this  shoulder  results 
from  an  avoided  crossing  with  a  state  with  Rydberg  character. 

The  first-derivative  coupling  between  the  X  and  A  states  peaks  at 
8.5  a^,  and  the  splitting  of  the  adiabatic  energy  values  is  a  minimum  at 
R  =  8.9  a^.  The  position  and  value  of  the  minimum  of  the  adiabatic  splitting 
is  compared  quantitatively  to  previous  calculations  in  Table  V.  The  table 
shows,  as  is  also  clear  from  Figure  1,  that  the  present  calculations  slightly 
overestimate  the  X-A  splitting.  However  the  present  calculations  are  more 
accurate  than  all  previous  calculations.  There  are  no  experimental  results 
available  for  the  other  avoided  crossings.  There  is  still  quite  a  bit  of 
uncertainty  about  the  splitting  at  the  A-C  avoided  crossing,  but  the  two 
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calculated  values  for  Che  minimum  splitting  of  the  X  and  C  curves  are  in 
reasonably  good  agreement. 

The  shape  of  the  present  diabatic  couplings  as  functions  of  R  are  very 
similar  to  those  calculated  previously  for  KH  (see  Figure  2  of  reference  7). 
The  adiabatic  couplings  also  show  a  reasonable  similarity  in  shape  to  those 
calculated  previously  (see  Figure  17  of  reference  7). 

E.  Semiclassical  approximations.  Although  F^2(R)  and  |e^(R)  - 
show  the  behavior  associated  with  an  avoided  crossing  at  R  =  8-9  a^.  Figure  3 
shows  that  the  diabatic  potential  curves  do  not  show  such  a  crossing.  Host 
previous  workers  have  treated  the  X-A  inelastic  transition  in  terms  of  such 
a  hypothetical  diabatic  crossing,  although  there  was  already  some  indication 
in  references  3  and  7  that  the  usual  diabatic  pictures  for  the  curve  crossings 
in  NaH  and  KH  might  be  inadequate.  Nevertheless  it  is  interesting  to  briefly 
compare  the  present  inelastic  transition  probabilities  to  those  obtained 
by  the  semiclassical  Landau-Zener  formula  for  estimating  transition  proba¬ 
bilities  at  adiabatic  avoided  crossings  resulting  from  diabatic  crossings. 
According  to  this  formula,  if  the  crossing  occurs  at  R  =  R^.  the  inelastic 
s-wave  transition  probability  is  ’ 


where 


V-  =  2p(l  -  p) 


-w 

p  =  e 


2rr 

=  fiv 


<Hu>2 


d(Hll  "  H22)/dR 


(51) 


(52) 

(53) 


and  v  is  the  local  radial  speed  at  the  diabatic  crossinn,  i.e.,  for  an  s  wave. 
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v  =  { 2 (E  -  H^1(Rx) l/p}4 


In  Che  two-state  approximation,  assuming  orthogonal  diabatic  states, 


H12(V  =  ’'5t£2(V  '  E1(V] 


(54) 


(55) 


5  6  8 

It  is  customary,  ’  ’  for  the  X-A  transition  in  alkali  hydrides,  to  neglect 
and  dH^/dR  at  R  =  and  to  approximate  by  (in  atomic  units) 


h22(e>  s  ip  •  “  -  R  - 


(56) 


where  IP  is  the  ionization  potential  of  M,  EA  is  the  electron  affinity  of  H, 
and  a  is  the  sum  of  the  polarizabilities  of  M+  and  H  ,  Putting  all  these 
approximations  together  and  noting  that  w  >>  1  yields  (in  atomic  units) 

,2 


P^  =  2  exp  [- 


tt|  Ae(Rx) 


(SE/p)^2  +  2aRxJ) 


(57) 


We  use  u  =  1790.83  a.u.  and  a  =  218  a.u.,  and  we  take  R^  and  !e(Rx)  from 
Table  V. 

In  the  more  complete  Landau-Zener-Stueckelberg  theory,  the  inelastic 
36,38-40 


probability  becomes' 


P12  =  2pLZ  sit,2(T  + 


(58) 


where  t  is  the  difference  in  action  integrals  for  the  two  adiabatic  curves 
and  6  is  a  phase  correction: 


Rv 


t  =  (2p/h“)  ‘\j  [E  -  e  (R) ]  'dr  -  j  it  -  t-iR) j  'dr;  U 

Rtpl  Rtp2 

where  R^^  and  R  ^  are  the  classical  turning  points  in  adiabatic  states  1 

LZ 

and  2  respectively.  Equation  (58)  shows  that  2P  should  be  compared  to  the 


upper  envelope  of  Che  oscillating  inelastic  transition  probability.  This 
comparison  is  shown  in  Table  VI  for  those  energies  at  which  the  close  coupling 
probabilities  are  maxima  (compare  Figure  5).  The  comparison,  however,  shows 
that  the  Landau-Zener  formula  overestimates  the  inelastic  transition  prob¬ 
abilities  as  compared  to  the  results  obtained  from  the  close  coupling  cal¬ 
culations.  The  overestimates  in  Table  VI  range  from  a  factor  of  36  to  a 
factor  of  60.  Thus  this  simple  high-energy  model  does  not  yield  accurate 
transition  probabilities  at  the  low  collision  energies  of  the  present  study. 

Although  the  failings  of  the  Landau-Zener  theory  are  well  known,  even 

41  42  43 

at  higher  energies,  both  Faist  and  Levine  and  Andresen  e£  al_.  found 

that  it  works  very  well  for  the  ionic-covalent  crossing  in  alkali-halogen 

43 

collisions,  even  near  threshold.  In  fact  Andresen  et_  al.  concluded  that 
"a  similar  agreement  is  expected  for  all  other  systems  which  are  dominated 
by  the  interaction  of  a  covalent  and  an  ionic  channel."  By  this  argument 
one  would  expect  that  it  would  be  accurate  for  the  ionic-covalent  inter¬ 
action  leading  to  In  previous  work,*7  however,  it  has  been  pointed 

out  that  the  strong  interaction  is  not  well  localized  in  this  case  and  that 
this  would  lead  to  a  breakdown  of  the  Landau-Zener  method.  Furthermore, 
the  fact  that  the  diabatic  curves  do  not  cross  is  an  indication  that  the 
transition  is  not  an  isolated  curve  crossing.  A  more  quantitative  measure 
of  whether  this  transition  should  be  treated  as  a  curve  crossing  will  be 
given  below. 

The  other  inelastic  transitions  are  more  complicated.  The  second  and 
third  diabatic  states  cross  twice,  and  the  first  and  third  diabatic  states 
do  not  cross.  Landau-Zener-tvpe  isolated-avoided-crossing  treatments  are 
not  appropriate  for  either  of  these  transitions.  In  the  absence  of  simple 
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models  Che  full  close  coupling  treatments  reported  above  are  the  best  way  to 
estimate  the  transition  probabilities. 

The  additional  factor  in  equation  (58)  results  from  the  interference 
of  the  two  possible  trajectories  leading  to  the  same  inelastic  collision. 
This  effect  accounts  for  the  existence  of  the  oscillations  in  Figure  5. 
Equation  (53)  results  from  a  high-energy  approximation,  retaining  only  the 

<  44 

leading  term  in  1/v  in  the  semiclassical  phase  integral.  With  the  same 
approximation  the  difference  in  action  integrals  becomes 


T 


R 


X 


(1/tiv)  / 
R 


Ae(R)  dr 


(60) 


tp 

where  the  classical  turning  points  are  the  same  to  this  order.  Equation  (60) 
predicts  that  at  high  energy  the  oscillations  should  be  evenly  spaced  in  1/v. 

At  the  low  energies  of  the  present  study,  the  initial  and  final  speeds  are 
appreciably  different,  and  they  differ  significantly  from  the  local  radial 
speeds  at  small  r.  Thus  it  is  not  possible  to  approximate  all  these  speeds 
by  the  same  v  and  equation  (60)  is  inapplicable.  Nevertheless  the  regularity 
in  the  oscillations  is  evident  in  Figure  5. 

We  have  already  mentioned  that  the  couplings  between  the  adiabatic  curves 
cannot  be  treated  as  simple  avoided  crossings.  To  put  this  on  a  more  quantitative 
basis  we  will  discuss  the  classification  of  the  nonadiabatic  couplings.  One 
may  distinguish  two  kinds  of  strong  interaction  in  terms  of  the  following 
two-state  representation 


<pf  (x,R)  =  <f>?(x)  cos  9(R)  +  4>,(x)  sinJ(R)  (61) 

^(x.R)  =  -t'J(x)  sin  3(R)  +  $*( x)cose(R )  (62) 


T 
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Using  (19),  (37),  (61),  and  (62),  one  finds 


fa  -  ii 

f12(R)  dR 


(63) 


One  can  distinguish  two  prototype  cases.  In  one  case,  6  changes  from  0  at 
R  =  ®  to  tt/4  (=0.785)  at  small  R.  This  occurs  for  a  symmetric  resonance 
interaction  like  H+  +  H.  In  another  case  9  changes  from  0  at  R  =  “  to 
■tt/2  (=1.57)  at  small  R.  This  occurs  for  a  diabatic  curve  crossing.  In 
general  one  may  associate  a  A6_  with  each  peak  in  f ^  (R)  by  integrating 
over  the  range  of  strong  interaction: 

R, 


A6 . .  =  /  fa.(R)  dR 
ij  p 

R1 


(64) 


For  actual  cases,  A9_  may  come  out  somewhere  between  0.785  and  1.57  for 
strong  interactions  and  may  be  less  for  weak  ones.  We  applied  this  model 
to  the  present  case  and  the  results  are  given  in  Table  VII.  These  results 
show  that  49^2  is  noC  close  to  1.57.  Thus  it  is  an  oversimplification  to 
describe  the  1-2  ionic-covalent  interaction  as  a  simple  curve  crossing, 
and  this  helps  to  explain  the  failure  of  the  Landau-Zener  method.  The 
adiabats  are  not  merely  "switching"  from  one  diabat  to  the  other;  rather 
there  is  an  appreciable  "mixing"  contribution.  Even  the  2-3  interaction, 
where  the  diabats  cross,  is  not  a  pure  curve  crossing.  Table  VII  shows 
that  the  1-3  interaction  is  weak. 

In  order  to  calculate  the  A0„  values  in  Table  VII,  we  had  to  separate 
the  overlapping  contributions  from  two  different  interaction  regions  in 
the  vicinity  of  the  sign  change  in  fa  (R) .  We  did  this  by  the  simplest 
method,  i.e.v  we  integrated  from  the  peak  of  fa^ (R)  to  R  =  »  and  multiplied 
by  2.  There  is  an  even  simpler  way  to  classify  the  nonadiabatic  interactions 


Chat  also  avoids  Che  question  of  overlapping  wings  of  the  peaks  in  f^.(R). 

Under  simple  assumptions  in  (R) ,  one  can  show  on  a  two-state  model  that 

ft j (R)  for  the  symmetric  resonance  has  the  form  (S/4)  sech  [S(R  -  R^) ]  and 

22 

for  the  curve  cross  cases  it  is  a  Lorentzian.  These  different  shapes 
may  be  characterized  by  defining  the  unitless  interaction  parameter: 

=  [peak  value  of  f^(R)]  *  FWHM 

where  FWHM  is  the  full  width  at  half  maximum  of  the  peak.  This  leads  to 

Q. .  =  0.66  for  the  svmmetric  resonance  case  and  0,.  =  1.00  for  the  curve 
ij  "ij 

crossing  cases.  We  also  calculated  for  the  biggest  peak  in  each 

f_(R),  and  the  results  are  in  Table  VII.  This  confirms  that  the  1-2 

and  2-3  interactions  involve  considerable  mixing  (as  in  symmetric  resonance) 


as  opposed  to  pure  curve  crossing.  This  simple  method  for  characterizing 
nonadiabatic  interactions  should  be  useful  for  many  problems. 
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VI.  CONCLUDING  REMARKS 

This  paper  has  been  concerned  with  calculational  techniques  for  the 
calculation  of  electronically  inelastic  transition  probabilities  using 
ab  initio  potential  energy  curves  and  nonadiabatic  radial-first-derivative 
couplings  as  input.  We  have  demonstrated  that  very  small  transition  prob¬ 
abilities  can  be  calculated  with  high  precision  using  a  convenient  and 
accurate  method  based  on  R  matrix  propagation  and  an  exponential  sector 
transformation  matrix. 

The  emphasis  in  this  article  is  on  the  techniques  for  solving  the 
coupled-channel  equations  to  obtain  the  precise  values  of  the  cross  sections 
that  correspond  to  a  set  of  al)  initio  potential  energy  curves  and  nonadiabatic 
radial-first-derivative  couplings.  Several  other  considerations  enter  when 
we  try  to  obtain  accurate  cross  sections,  i.e.,  good  agreement  with  experi¬ 
ment  or  reliable  predictions.  Only  two  of  these  have  been  considered  in 
this  paper,  namely,  sensitivity  to  changes  in  the  potential  curves  and  to 
the  treatment  of  the  second-derivative  nonadiabatic  coupling  terms.  Some 
other  factors  that  must  be  considered  in  future  work  to  obtain  reliable 

results  for  the  K  +  H  system  are  sensitivity  to  change  of  electronic  origin 

3  22  A  5 

in  the  radial-first-derivative  coupling  terms  '  ’  and  the  role  of  angular- 

22  A  6 

derivative  coupling  terms.  ’  The  latter  terms  couple  the  I  states  to  the 

tl  -  Mates ,  but  we  have  included  only  Z—Z  coupling  in  the  present  study. 

2  3  47  48 

Another  question,  also  discussed  elsewhere  ’  *  ’  and  also  deferred  to 

future  work  for  detailed  numerical  study,  is  the  question  of  whether  accurate 
low-energy  cross  sections  can  be  calculated  from  ab  initio  molecular-frame 
input  data  without  either  including  plane-wave  factors  or  transforming  at 
large  R  to  a  laboratory-frame  diagonal  representation  in  which  the  coupling 


d 


vanishes  at  infinity  for  any  consistent  choice  of  origin  for  the  electronic 
coordinates.  If  such  a  transformation  is  required,  the  present  method 
can  still  be  used  to  integrate  out  to  the  large-R  transformation  distance. 
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APPENDIX 

In  this  appendix  we  discuss  the  fact  that  using  (18)  with  an  incomplete 

basis  is  equivalent  to  including  the  full  effect  of  the  first  derivative 
£ 

coupling  matrix  F  (R) ,  but  including  only  part  of  the  second  derivative 

Vi 

2 

coupling  matrix  G  (R) . 

First  note  that  U(R)  defined  by  (30)  and  M(R)  defined  by  (18)  obey  the 

same  differential  equation.  Note  that  the  sign  difference  arises  in  taking 

the  transpose  of  (30)  since  fa(R)  is  skew  symmetric.  It  has  been  shown  else- 
3 

where  that  the  matrices  M(R)  can  be  used  to  construct  the  transformation 
U(R)  from  the  adiabatic  basis  to  the  diabatic  one.  Therefore,  the  two 
problems  are  equivalent  in  the  limit  of  small  sector  size  and  since  the 
transformation  method,  (27)-(34),  includes  the  full  effect  of  the  Fa(R) 
matrix,  so  will  the  direct  integrationscheme  of  eqs.  (7)-(26). 

By  using  the  M  matrices  to  generate  the  transformation  to  the  diabatic 

basis  \^(R),  (7)  can  be  transformed  to 

v 

2 

xd(R;i)  =  D(R;i)  xV;i)  (A-l) 

dR“  ^  % 

where 

D(R;i)  =  U(R^)  ^2(R;i)  UT(Rj)  (A-2) 

and 

Xd(R;i)  =  U(Ra)  ya(R;i)  (A-3) 

r\j 

In  enforcing  the  sector  matching  conditions  (10)  and  (11),  the  matrix  T(i-l;i) 
is  now  replaced  by  the  identity  matrix  since  Fd(R)  is  identical  zero.  There¬ 
fore,  solution  of  (A-l)-(A-3)  is  the  same  as  the  solution  to  (31)  but  neglecting 
Gd (R) .  Since  G^(R)  vanishes  in  the  limit  of  a  complete  set,  this  neglect 


is  completely  justified  in  the  converged  limit.  Consider,  however  the  case 
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in  which  an  incomplete  basis  of  the  N  lowest-energy  adiabatic  states  is 
included  in  the  coupled-channel  equations.  Then  the  terms  neglected  in 
assuming  (R)  obey  (35)  are  given  by 


VR) 


AB  k=N+l 


l  <♦!! 


3  I  o  \  3  a  |  8  i  3 
<6  .  rnrl  4>,  ><(t> ,  !— !  ? .  > 
a  3R'*k  k 1  aR  '  v l 


(A-4) 


Equation  (35)  will  be  a  valid  approximation  for  incomplete  bases  in  the  case 


f(R)  +  Ga(R)i  »  |  ACR) 


(A-5) 


Using  a  Hellman-Feynman  type  theorem  it  can  be  sliown  that 


i  3H  1 

?a|^!pf>  =  - i -  a 

1  jR'  k  ea(R)  -  ea(R)  1  L'R  k 
k  1 


(A-6) 


and  therefore  the  conditions  for  the  validity  of  (35)  are 

3H  i  „ 
<^j_ei|ia>2 

i  ‘  k 


>2 

£i(R)  +  GuCR)i  '  > 


2wAB  'k-N+1  [-t  J(R)  -  ta(R)  ]~  : 


(A- 7) 


and 


|G®, (R) !  » 


ij 


AB 


l 


3H  i  3H  i 

a |  e1! ,a  a i  ex <  a 

<Pil-3riV<7k:  — !^> 


k=N+l  [ea(R)  -  ta(R)][-;a(R)  -  eJ(R) 


(A-S) 


Thus  equation  (35)  is  a  good  approximation  if  the  omitted  states  lie  high 


enough  in  energy. 


TABLE  I.  Relative  computing  times  required  to  obtain  0.5 %  accuracy 


Method 

Computing  time 

D. 

R  matrix  propagation 

1.0 

C. 

Magnus,  diabatic,  analytic  e::t 

onential 

**  r\ 
t-  •  yj 

B. 

Magnus,  diabatic,  exponential 

by  power  series 

4 

A. 

Magnus,  adiabatic 

7 

G. 

Runge-Kutta,  diabatic 

16 

E . 

Runge-Kutta,  adiabatic 

16 

F. 

Adams-Moulton,  adiabatic 

>30a 

H. 

Adams-Moulton,  diabatic 

>40a 

In  these  two  cases  the  results  are  still  accurate  to  only  10-50%  for 
the  computing  times  listed.  Since  these  two  methods  were  found  to  be 
so  inefficient,  we  did  not  continue  to  decrease  the  stepsize  parameter 
to  obtain  0.5%  accuracy. 


i 
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unclassified 


TABLE  II.  Detailed  comparison  of  computational  efficiencies  of  methods 


C  and  D. 


Number  of  sectors 
or  steps 

Computing 
time  a 

%  error 

Method  D 

302 

0.85 

1. 5646 (—6) C 

0.68 

358 

1.01 

1.5617 (-6) 

0.49 

452 

1.23 

1 . 5564 (-6) 

0.37 

486 

1.33 

1.5554 (-6) 

0.09 

520 

1.43 

1 . 5550(-6) 

0.06 

621 

1.68 

1 . 5543 (-6) 

0.02 

659 

1.73 

1.5541(-6) 

0.01 

926 

2.56 

1 . 5539 (— 6 ) 

-0.01 

2344 

6.25 

1 . 5539 (-6) 

-0.01 

4671 

11.82 

1.5540(-6) 

0.00 

Method  C 

453 

1.52 

1.5741 (-6) 

1.29 

647 

2.02 

1 . 5628 (-6) 

0.57 

653 

2.03 

1.5615(-6) 

0.48 

680 

2.09 

1 . 5612 (-6) 

0.46 

1359 

3.55 

1.5558(-6) 

0.12 

1477 

4.66 

1 . 5555 (-6) 

0.10 

4527 

11.93 

1 . 5541 (-6) 

0.01 

13531 

35.76 

1. 5540(-6) 

0.01 

Same  scale  as  Tablj  I. 


^For  2H  potential  curves,  S  =  0.06  E,  ,  i  = 

n 

CNumbers  in  parentheses  are  powers  of  ten. 


0. 


TABLE  III.  Comparison  of  s-wave  inelastic  transition  probabilities  calculated  with  two  different  equations 


Numbers  in  parentheses  are  powers  of  ten. 


TABLE  IV.  Cross  sections  for  excitation  and  de-excitation  processes 
in  K  +  H  collisions.3 


E<Eh> 

i  b 
max 

°12  (a0) 

°21  (a0>  C’e 

°13  (a0) 

,  2%e 
°23  (a0} 

21  potential  curves 

0.06 

37 

8.40(-6)d 

2. 08 (-4) 

0.10 

165 

1. 20(-3) 

9.80(-4) 

2. 75(-6) 

0.41 

2H  potential  curves 

0.06 

36 

1.60(-5) 

3.96(-4) 

0.10 

>145 

9. 24 (-4) 

7.55 (—4) 

1. 73(-6) 

0.4 

Cross 

sections  are  accurate  to 

i  about  2%  with 

respect  to  variations 

of  integration  parameters. 

3Number  of  partial  waves  necessary  to  converge  inelastic  cross  sections 
to  5  significant  figures. 

-  2  2 

uo..  =  (k.d . /k.d  .  )o_,  .  where  lik.  is  asymptotic  momentum  in  channel  i 
J1  i  i  J  y  ij  1 

and  d^  is  the  degeneracy  of  state  i. 

^Numbers  in  parentheses  are  powers  of  ten. 
eIncludes  factor  of  1/3  for  P-state  degeneracy. 


TABLE  V.  Predicted  positions  R  and  values  Ae  of  the  large-R  minima  in  the  splitting  of  the  adiabatic 


TABLE  VI.  Maxima  of  s-wave  inelastic  transition  probabilities 


for  X-*-A  transition. 


E 

(eV) 

P12 

close  coupling 

2PLZ 

Method 

21 

1.91 

1.6(-5)3 

8.4(-4) 

2.22 

3.3(-5) 

1  •  5 (-3) 

2.96 

8 • 9 (-5) 

4 .4  (-3) 

Method 

2H 

2.00 

2 . 8 (-5) 

1.0(-3) 

2.34 

3.2(-5) 

1.9(-3) 

3.11 

1.4 (-4) 

5 . 2 (-3) 

Numbers  in  parentheses  are  powers  of  ten. 
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Figure  Captions 


Fig.  1.  Adiabatic  potential  energy  curves  as  a  function  of  internuclear 

distance  for  the  three  lowest  '*T+  states  of  KH.  The  curves  are  the 

results  of  the  ab  initio  pseudopotential  calculations  as  obtained  by 

method  21  for  the  solid  curve  and  method  2H  for  the  dashed  curve. 

The  two  asymptotic  values  for  the  C  state  differ  because  in  fitting 

2 

the  21  potential  curve  the  experimental  5  S  excitation  energy  of 

2.61  eV  was  used,  whereas  for  the  2H  potential  curve  the  numerically 

computed  value  of  2.55  eV  was  used.  This  difference  has  a  negligible 

effect  on  the  scattering  calculations.  The  potential  curves  for  the 

2 

B  state  both  dissociate  to  the  correct  atomic  4  P  excitation  energy 
of  1.61  eV.  The  points  are  the  experimentally  determined  RKR  values 
for  the  X  and  A  potential  curves  (from  references  7  and  10). 

Fig.  2.  First  derivative  coupling  terms,  as  defined  by  equations  (19)  and 
(37),  for  the  three  lowest  adiabatic  states  of  KH  as  functions 
of  internuclear  distance.  States  1,  2,  and  3  correspond  to  the  X, 

A,  and  C  states  of  figure  1.  (Since  the  difference  between  methods 
21  and  2H  involves  only  the  treatment  of  the  core,  the  derivative 
couplings  are  the  same  for  both  methods.) 

Fig.  3.  The  diagonal  matrix  elements  H^(R)  for  the  ^E+  states  of  KH  in 

the  P-diabatic  basis  as  functions  of  internuclear  distance.  The  solid 


and  dashed  curves  are  the  results  of  transforming  the  corresponding 
solid  and  dashed  adiabatic  curves  shown  in  figure  1. 


Fig.  4.  Same  as  figure  3  except  for  the  off-diagonal  elements  H.,  (R). 

J  K 

Fig*  5.  Transition  probabilities  P.,(E)  as  functions  of  the  reciprocal 

3 

of  the  total  energy  for  ?.  =  0.  The  zero  of  energy  is  the  asymptote 


of  che  X  state.  The  top  plot  is  for  excitation  of  K  from  the  4  S 


to  the  4  P  state,  the  lower  left  plot  is  for  excitation  of  K  from 
2  2 

the  4  S  to  the  5  S  state,  and  the  lower  right  part  is  for  excitation 
2  2 

of  K  from  the  4  P  to  the  5  S  state.  In  the  top  plot  the  left  and 

right  arrows  along  the  abscissa  indicate  the  energetic  thresholds 

2  2 

for  excitation  to  the  5  S  and  4  P  state,  respectively.  The  arrows 

in  the  lowest  two  plots  are  the  threshold  for  the  excitation  to  the 
2 

5  S  state.  Note  that  the  energy  scale  for  the  upper  and  leer  left 
plots  coincide  and  are  aligned  vertically  with  each  other,  whereas 
the  energy  scale  for  the  lower  right  plot  is  expanded.  In  each  plot 
the  solid  and  dashed  curves  are  the  probabilities  computed  using 
the  adiabatic  potential  curves  and  nonadiabatic  couplings  obtained 
by  methods  21  and  2H,  respectively. 
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